Mechanisms of Qu's Formula 3 for Improving Endometrial Receptivity in PCOS: A Combination of Network Pharmacology and Molecular Dynamics Simulation

Summary

This protocol aims to explore the mechanism by which QUF3 improves endometrial receptivity in PCOS using network pharmacology, docking, and MD simulation.

Abstract

Impaired endometrial receptivity is a major cause of infertility in polycystic ovary syndrome (PCOS). Qu's Formula 3 (QUF3) is a Chinese herbal medicine formula clinically used to improve endometrial receptivity in patients with PCOS. This study employed network pharmacology, molecular docking, and molecular dynamics (MD) simulations to investigate the underlying mechanisms. Active constituents of QUF3 were identified using the TCMSP database, and potential targets related to endometrial receptivity and PCOS were retrieved from DrugBank and other resources. A compound-target interaction network and a protein-protein interaction (PPI) network were constructed via Cytoscape to identify key targets. Core targets were subjected to GO and KEGG enrichment analyses. Molecular docking, MD simulations, principal component analysis (PCA), free energy landscape (FEL), and dynamic cross-correlation matrix (DCCM) were used to evaluate binding interactions. From 91 active ingredients and 294 potential drug targets, 60 disease-related targets were identified. Luteolin and sesamin were among the key pharmacodynamic components. Ten core targets were identified: AKT1, EGFR, TNF, TP53, IL6, BCL2, ESR1, IL1B, STAT3, and MMP9. KEGG enrichment revealed 132 signaling pathways, and GO analysis identified 678 entries. MD simulations indicated that the binding between the top five active constituents and their respective targets was stable. PCA, FEL, and DCCM further demonstrated high thermodynamic stability and structural rigidity of these complexes. In conclusion, QUF3 improves endometrial receptivity in PCOS through multicomponent, multitarget, and multipathway interactions. This study provides a theoretical basis, from a computational simulation perspective, for the development of targeted Chinese herbal medicine therapies for PCOS-related infertility, and may have positive implications for improving pregnancy outcomes in women with PCOS in the future.

Introduction

Polycystic ovary syndrome (PCOS) is one of the most common endocrine disorders in women¹. Its main characteristics include chronic persistent ovulatory dysfunction, clinical or biochemical hyperandrogenemia, and polycystic ovarian morphology, with clinical manifestations such as amenorrhea, infertility, hirsutism, acne, and obesity. In recent years, advances in ovulation induction and assisted reproductive technology (ART) have significantly improved pregnancy rates in PCOS patients; however, pregnancy outcomes remain poorer than those in healthy women^2,3. Endometrial receptivity is a key factor influencing ART outcomes⁴. Endometrial receptivity refers to the ability of the endometrium to accept an embryo at a specific physiological stage, which plays a decisive role in embryo implantation and subsequent development^5,6. Pathological conditions commonly observed in PCOS, such as obesity, insulin resistance, and chronic inflammation, can impair endometrial receptivity, hinder fertilized egg implantation and embryo development, and ultimately lead to infertility and miscarriage^7,8. Nevertheless, the regulatory mechanisms underlying endometrial receptivity in PCOS remain largely unclear, and the interaction network among various influencing factors has yet to be elucidated. Qu's Formula 3 (QUF3), a Chinese herbal medicine formula clinically used to improve endometrial receptivity in patients with PCOS, is supported by solid classical theories of traditional Chinese medicine (TCM) in gynecology. The formula consists of Paeoniae Radix Alba (Baishao), Rehmanniae Radix Praeparata (Shudihuang), Cuscutae Semen (Tusizi), Fructus Ligustri Lucidi (Nvzhenzi), Herba Taxilli (Sangjisheng), Radix Salviae (Danshen), and Cornus Officinalis Sieb. et Zucc. (Shanzhuyu) in a ratio of 2:4:3:4:3:2:2, prepared by ethanol reflux extraction. QUF3 exerts synergistic effects of multiple herbs, tonifying the kidney, enriching essence, nourishing blood, and promoting blood circulation⁹. Its clinical efficacy in improving endometrial receptivity in PCOS patients has been well validated, but its core pharmacological mechanisms remain to be further investigated.

Network pharmacology is a systems biology-based analytical approach that constructs multi-level “herb-component-target-disease” networks to reveal, from a holistic perspective, the intervention effects of drugs on disease networks. Unlike the traditional “single drug-single target” model, network pharmacology can comprehensively capture the synergistic characteristics of "multi-component, multi-target, and multi-pathway" actions of TCM formulas, making it particularly suitable for investigating the complex mechanisms of TCM. Molecular dynamics (MD) simulation, in contrast, is a computational technique based on Newtonian mechanics that simulates, at the atomic level, the conformational changes and dynamic behavior of protein-ligand complexes over time. Compared with static molecular docking, MD simulation can assess dynamic parameters, such as binding stability, free energy landscapes, and structural rigidity, thereby providing a more realistic evaluation of interactions between drug molecules and target proteins. Network pharmacology is suitable for globally screening potential targets and pathways of TCM formulas, while MD simulation can dynamically evaluate the binding stability between key components and target proteins. Combining these two approaches enables a more reliable analysis of the multi-component, multi-target synergistic mechanisms of TCM formulas. Based on this rationale, this study employed a combined network pharmacology and MD simulation approach to systematically explore the potential pharmacological mechanisms underlying QUF3's improvement of endometrial receptivity in PCOS patients.

This combined computational workflow is most suitable for the initial exploration of the multi-component, multi-target mechanisms of Chinese herbal formulas or natural products, particularly for complex systems in which the active ingredients and targets are largely unknown. Users should have basic bioinformatics knowledge, including database searching and network analysis, as well as the ability to perform fundamental computational biology analyses such as molecular docking and molecular dynamics simulations using common software packages. A Linux-based computing environment with at least 16 GB of RAM is recommended, along with access to relevant public databases such as the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP), PubChem, and UniProt. All software used in this workflow (e.g., Cytoscape, AutoDockTools, GROMACS) is open source or free for academic use, and readers can obtain it from the information provided in the Table of Materials.

Protocol

Screening of active constituents and potential targets of QUF3
The active chemical components of Baishao, Shudihuang, Tusizi, Nvzhenzi, Sangjisheng, Danshen, and Shanzhuyu in the medication QUF3 were identified using the TCMSP on April 1, 2025¹⁰. The active ingredients were preliminarily screened based on the conditions of oral bioavailability (OB) ≥30% and drug-likeness (DL) ≥ 0.18¹¹. These thresholds were chosen as standard empirical criteria to ensure sufficient absorption and drug-like properties. The canonical SMILES sequence of each compound was searched in the PubChem database, and the resulting sequence was used to predict the compound targets in the Swiss Target Prediction online database, followed by the collection of target protein gene names and UniProt IDs from the UniProt database^12,13,14. As an intermediate checkpoint, the number of unique targets obtained after deduplication was recorded.

Identification of action targets for endometrial receptivity in PCOS
The keywords “polycystic ovarian syndrome” and “endometrial receptivity” were used as search terms in the Online Mendelian Inheritance in Man (OMIM), DrugBank, Therapeutic Target Database (TTD), and GeneCards databases to identify potential targets^15,16,17,18. For GeneCards, only targets with a relevance score >5 were retained; for the other databases, all retrieved targets were included. Following merging and removal of duplicate entries, the number of unique disease-associated targets was recorded as an intermediate output.

Network construction of components and targets
The intersection targets were generated by entering the possible targets for QUF3, PCOS, and endometrial receptivity into an online Venn diagram tool. These intersection targets were considered as the key targets through which QUF3 exerts its efficacy in improving endometrial receptivity in PCOS. Subsequently, a “herb-constituent-target” (H-C-T) network was created by matching these intersection targets with the active constituent-target data of QUF3 using Cytoscape¹⁹. The network was visualized using the “Network Analyzer” tool.

Protein-protein interaction (PPI) network construction
Potential target genes of QUF3 for improving endometrial receptivity in PCOS were submitted as input to the STRING data platform²⁰. The species was set to "Homo sapiens," and the confidence threshold was set to “medium confidence (0.4)”, yielding a TSV file of protein-protein interaction data. This file was then imported into network analysis software to construct and analyze the PPI network²¹. Node importance in the network is indicated by the “Degree” value: a higher Degree value indicates a more significant node. Subsequently, the primary active ingredients were sorted by Degree, and the top 10 core targets were selected based on Degree ranking.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes pathway enrichment (KEGG) analysis
GO and KEGG pathway enrichment analyses were conducted using the DAVID database to explore the biological pathways and mechanisms by which QUF3 improves endometrial receptivity in PCOS²². The species was set to “Homo sapiens”, and the significance threshold was set at P < 0.05. The GO analysis covered biological process (BP), cellular component (CC), and molecular function (MF)²³. The results were sorted by the number of enriched targets. Pathways related to human metabolism and diseases were then excluded, and the top 20 KEGG pathways were selected for visualization. A “component-disease-KEGG pathway” multidimensional network was subsequently constructed using network analysis software.

Molecular docking analysis
The top 10 key active constituents and the top 10 core targets were selected for molecular docking analysis using AutoDockTools and PyMOL software^24,25. PDB files of the key targets were downloaded from a protein structure database²⁶, and SDF files of the main active components were obtained from a chemical database²⁷. The protein molecules were then preprocessed using molecular docking software, which included desalting, hydrogenation, and charge calculation. Semi-flexible docking was performed using a grid box covering the entire protein-binding site. The docking process was considered successful when the binding energy was negative, and the root-mean-square deviation (RMSD) of the top pose was ≤ 2.0 Å. After docking, the results for the main active components of QUF3 with the core targets were visualized using molecular visualization software.

Molecular dynamics (MD) simulations
MD simulations of the protein–ligand complex were performed using GROMACS software to explore the interaction between receptors and ligands^28,29. The amber99sb-ildn force field and the general Amber force field (GAFF) were used to generate the parameters and topologies of proteins and ligands, respectively. The size of the simulation box was set such that the distance between each protein atom and the box edge was greater than 1.0 nm. The box was filled with an explicit solvent using the simple point charge model (SPC216 water molecules), and water molecules were replaced with Na⁺ and Cl⁻ counterions to make the simulation system electrically neutral. The entire system was optimized by the steepest descent method to reduce unreasonable contacts or atom overlaps. Sufficient pre-equilibration of the simulation system was achieved by performing NVT and NPT ensembles for 100 ps each at 300 K and 1 bar, respectively. Subsequently, a 100 ns MD simulation was performed under periodic boundary conditions, with temperature (300 K) and pressure (1 bar) controlled by the V-rescale and Parrinello-Rahman methods, respectively³⁰. The Newtonian equations of motion were solved using the leapfrog integrator with a time step of 2 fs. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method with a Fourier spacing of 0.16 nm, and all bond lengths were constrained using the LINCS method. The VMD software was used to display, analyze, and animate trajectories³¹. The binding free energy of each compound was calculated using a binding free-energy calculation script.

Principal component analysis (PCA), free energy landscape (FEL), and dynamical cross-correlation matrix (DCCM) analysis
PCA was conducted using MD simulation software to analyze the dominant conformational motions of apo and ligand-bound complexes. Covariance matrices of Cα atomic fluctuations were generated, from which eigenvectors (collective motions) and eigenvalues (motion magnitude) were obtained³². Key dynamic regions and structural flexibility associated with ligand binding were identified through this analysis. Molecular trajectories were projected onto the first two principal components (PC1 and PC2) using the gmx anaeig tool, thereby capturing essential motions³³. Only the equilibrated final 100 ns (RMSD-stable phase) was used for PCA and FEL calculations. FEL was derived from the PCA data using gmx sham to compute the Gibbs free energy³⁴, revealing stable conformations and energy minima that reflect conformational diversity and state transitions. Dynamical cross-correlation matrix (DCCM) analysis of Cα atoms was performed to evaluate structural dynamics and inter‑residue motion coupling³⁵. This analysis elucidated the stability, intermolecular interactions, and collective motions of the protein-ligand complexes. The cross-correlation coefficient Cij between residues i and j was calculated as:

In the formula, the letters i and j represent the Cα atoms of two residues. When the value of Cij is positive and greater than zero; the two atoms moving in the same direction are considered to be in correlated motion. When the value of Cij is negative and less than zero, the two atoms moving in opposite directions are considered to be in anti-correlated motion. A Cij value of zero indicates no relationship between the motions of the two atoms. Furthermore, Cij values greater than 0.7 were interpreted as strongly correlated (positive) or strongly anti-correlated (negative) motions. To verify convergence, the RMSD and Rg profiles were examined. The stable profiles confirmed that the final 100 ns was sufficient for DCCM analysis. These motion patterns revealed how ligand binding modulates intramolecular communication and domain rearrangements. To ensure dynamic stability and convergence, only the final 100 ns of equilibrated trajectories were analyzed, and correlation maps were visualized to highlight regions with significant dynamic interplay.

To better explain the interaction energy between proteins and ligands, a binding free-energy calculation script was used to calculate the binding energy of all protein-ligand complexes at equilibrium. In the application of the binding free-energy calculation method³⁶, the total binding energy was decomposed into four independent parts: electrostatic interaction, van der Waals interaction, polar solvation, and non-polar solvation interaction. The non-polar solvation term is usually referred to as SASA.

Concluding procedural endpoint for target and pathway selection
At the end of the workflow, final validated targets and pathways were selected based on the following criteria: core targets with the top 10 Degree values in the PPI network, and key active constituents with the top 10 Degree values in the H-C-T network. For KEGG pathways, after excluding those related to human metabolism and diseases, the top 20 pathways were ranked and selected based on the number of enriched targets. MD simulation results were used to confirm stable binding between top constituents and core targets. These steps collectively defined the protocol's final output.

Results

Screening of active constituents and potential targets of QUF3
Using OB ≥ 30% and DL ≥ 0.18 as screening criteria, after searching the TCMSP database, QUF3 identified 126 chemical components. Among these, Paeoniae Radix Alba (Baishao) contained 13 active constituents, Rehmanniae Radix Praeparata (Shudihuang) contained 2, Cuscutae Semen (Tusizi) contained 13, Fructus Ligustri Lucidi (Nvzhenzi) contained 13, Herba Taxilli (Sangjisheng) contained 2, Radix Salviae (Danshen) contained 65, Cornus Officinalis Sieb. Et Zucc. (Shanzhuyu) contained 20 (as shown in Supplementary Table 1). The 126 chemical components were entered into the Swiss Target Prediction Online Database. After removing duplicates, 91 chemical components and 284 targets were obtained.

The targets of endometrial receptivity in PCOS
By searching the DrugBank, GeneCards, TTD, and OMIM databases, a total of 2847 disease targets were identified. Among them, there were 1800 disease targets for PCOS and 1047 targets related to endometrial receptivity. After removing duplicates, 348 disease targets were obtained.

Construction of H-C-T network
An intersection analysis of the disease targets and the potential target sites of QUF3 identified 60 common targets. Visual processing was carried out by constructing the Venny (Figure 1A). The intersection targets were then imported into Cytoscape to construct the H-C-T network. The network comprised 158 nodes and 326 edges, with 91 active drug constituents corresponding to 60 common targets (Figure 1B). Using the Network Analyzer and CytoNCA plugin, parameters such as connectivity, degree, and closeness were analyzed to identify active ingredients, including quercetin, luteolin, kaempferol, tanshinone IIA, beta-sitosterol, cryptotanshinone, 2-isopropyl-8-methylphenanthrene-3,4-dione, danshenxinkun d, sesamin, and isorhamnetin. The key active ingredients are shown in Table 1 and are predicted to play an important role in the QUF3 intervention process on endometrial receptivity in PCOS.

PPI network
The PPI network was constructed using 60 intersection targets (Figure 1C), yielding 60 nodes and 924 edges for subsequent analysis. Based on the degree of freedom, the core proteins were sorted out and 10 of them were selected, namely AKT serine/threonine kinase 1 (AKT1), epidermal growth factor receptor (EGFR), tumor necrosis factor (TNF), tumor protein 53 (TP53), Interleukin-6 (IL6), B-cell lymphoma-2 (BCL2), estrogen receptor 1 (ESR1), interleukin-1beta (IL1B), signal transducer and activator of transcription 3 (STAT3), and matrix metalloproteinase 9 (MMP9). These proteins were identified through computational network analysis and proposed as potential key targets for QUF3 to improve endometrial receptivity in PCOS.

GO and KEGG analyses
The GO enrichment analysis encompassed 678 GO terms, comprising 532 BP, 40 CC, and 106 MF. Figure 1D illustrates the top ten significantly enriched terms in the BP, MF, and CC categories. BP was mainly enriched in positive regulation of gene expression, transcription by RNA polymerase II, miRNA transcription, DNA-templated transcription, cell population proliferation, and the MAPK cascade. CC was mainly enriched in the protein-containing complex, extracellular region, extracellular space, platelet alpha granule lumen, and caveola. MF was mainly enriched in enzyme binding, identical protein binding, cytokine activity, and growth factor activity.

A total of 132 KEGG pathway enrichment items were identified. After excluding pathways related to human metabolism and diseases, the top 20 pathways were ranked by the count of enriched targets and selected for visualization (Figure 1E). These pathways converged on several classical signaling axes, including PI3K-Akt, MAPK, HIF (both HIF-1 and HIF-2), FoxO, JAK-STAT, and Ras/Rap1 signaling, as well as processes such as cellular senescence, focal adhesion, thyroid hormone and relaxin signaling, cytokine-cytokine receptor interaction, and immune-related pathways (IL-17, TNF, T cell receptor, and C-type lectin receptor signaling), along with osteoclast differentiation and prolactin signaling. Collectively, these pathways suggest that QUF3 may improve endometrial receptivity by modulating inflammation, metabolism, hormone responses, and cellular senescence. A “component-target-KEGG pathway” network was constructed in Cytoscape 3.10.3 (Figure 1F).

Molecular docking analysis
Using AutoDock, the active sites of the key target proteins were predicted (as shown in Supplementary Table 2). In this study, a binding energy ≤–5.0 kcal/mol was considered indicative of strong binding activity. Generally, the lower the binding energy value (in kcal/mol), the more stable the binding³⁷. As shown in Figure 1G, the active components Luteolin, Tanshinone IIa, Danshenxinkun d, and Sesamin in QUF3 had strong binding activity with the core targets, AKT1, MMP9, and BCL2. For example, the binding energy of luteolin-MMP9 was −10.3 kcal/mol, and that of sesamin-AKT1 was −10.6 kcal/mol.

Molecular docking visualization analysis was conducted on the top five active components and core targets using PyMOL. Sesamin formed a stable binding conformation with MMP9, with a binding energy of −11.9 kcal/mol (Figure 2A). Danshenxinkun D tightly bound to AKT1, showing a binding energy of −11.6 kcal/mol (Figure 2B). 2-isopropyl-8-methylphenanthrene-3,4-dione docked into AKT1 with a binding energy of −11.1 kcal/mol (Figure 2C). The Sesamin-AKT1 complex exhibited a binding energy of −10.6 kcal/mol, and hydrogen bond interactions were observed (Figure 2D). Luteolin-MMP9 had a binding energy of −10.3 kcal/mol, and this complex was stabilized by multiple non-covalent interactions acting cooperatively (Figure 2E). These visualization results further support the binding energy data. MD simulation was conducted to explore potential binding interactions. These computational results indicate that the active components may improve endometrial receptivity.

Molecular dynamics (MD) simulations
In MD simulations, RMSD analysis assesses structural stability: lower RMSD values indicate greater conformational stability, while higher values indicate increased flexibility. As shown in Figure 3A, all protein-ligand complexes reached equilibrium after 20 ns with stable RMSD profiles. Protein-ligand recombination SASA values did not vary substantially throughout all recombination simulations, as illustrated in Figure 3B, suggesting consistent protein-ligand binding. Radius of gyration (Rg), a key metric for evaluating structural compactness, reflects the spatial distribution of atoms relative to the molecular axis and assesses the stability of protein-ligand complexes. As shown in Figure 3C, Rg values remained stable throughout the MD simulation, which is consistent with SASA results and suggests robust protein-ligand binding. Protein-ligand interactions are stabilized by hydrogen bonding, which affects complex stability and binding affinity. As shown in Figure 3D, the hydrogen bond numbers of MMP9-Sesamin, AKT1-Sesamin, AKT1-Danshenxinkun D, MMP9-Luteolin, and AKT1-2-isopropyl-8-methylphenanthrene-3,4-dione complexes were 1, 3, 4, 5, and 2, respectively, suggesting that the protein and ligand interacted through hydrogen bonds.

Principal component analysis (PCA), free energy landscape (FEL), and dynamical cross-correlation matrix (DCCM) analysis
PCA is a commonly used unsupervised dimensionality reduction method that converts high-dimensional data into a low-dimensional space to observe protein stability and function³⁸. And the FEL delivers key insights into the mechanistic basis of protein folding by mapping conformational state transitions during native refolding³⁹. As shown in Figure 4A,B, the PCA and FEL suggest a stable combination of protein and protein-ligand complex. The DCCM analysis was performed for five protein-ligand complexes (Figure 4C), revealing distinct motion coupling patterns. The luteolin-MMP9 complex exhibited the strongest intradomain positive correlations (red patches), consistent with its high number of stable hydrogen bonds. For AKT1 complexes, 2-isopropyl-8-methylphenanthrene-3,4-dione induced the most extensive long-range anti-correlated motions (blue regions), suggesting robust allosteric regulation. Overall, the DCCM results suggested that ligand binding modulated protein dynamics in a compound-specific manner, with stronger correlations associated with higher binding stability.

The binding free energy analysis was performed to characterize the stability and binding affinity of the target protein–ligand complexes (as shown in Table 2). In the protein-ligand complex system, the free energy of protein-ligand binding for the MMP9-Sesamin, AKT1-Sesamin, AKT1-Danshenxinkun D, MMP9-Luteolin, and AKT1-2-isopropyl-8-methylphenanthrene-3,4-dione complexes was −184.853, −91.035,−79.089, −50.300, and −48.856 kJ/mol, respectively, indicating that the protein-ligand binding is stable and mainly involves van der Waals interactions.

DATA AVAILABILITY:
The raw data presented in this study are fully open and publicly available for anyone to download. All raw data have been deposited in the Science Data Bank and can be accessed via the DOI: 10.57760/sciencedb.36552 (https://doi.org/10.57760/sciencedb.36552).

Figure 1: The Network Pharmacology of QUF3 improves endometrial receptivity in PCOS. (A) The targets of compounds and disease. Figure legends: PCOS: polycystic ovarian syndrome; ER: endometrial receptivity. (B) The herb-constituent-target network diagram. Figure legends: The blue hexagon represents the target, the purple square represents the component, and the orange circle represents the drug. (C) The PPI analysis. The greater the degree value, the larger the circle, the redder the color; the smaller the degree value, the smaller the circle, the bluer the color. (D) The GO functional enrichment analysis. The abscissa represents GO functional categories, from left to right: biological processes, cell components, and molecular functions; the ordinate represents enrichment scores; the top 10 GO functional categories with the smallest P-values were selected. (E) The KEGG pathway enrichment analysis. The abscissa represents the Enrichment score, and the ordinate represents the pathway. The bubble area reflects the number of genes contained in the pathway. The larger the number of genes, the larger the bubble area. The color of the bubble indicates the size of the P-value; the smaller the P-value, the redder the color. (F) The compound-target-KEGG passway network diagram. The purple arrow represents the pathway, the blue hexagon represents the target, and the orange square represents the component. (G) Heat map of the docking between active ingredients and target proteins. Please click here to view a larger version of this figure.

Figure 2: Molecular docking visualization of the top five active constituents from QUF3 with their corresponding core targets. (A) Sesamin-MMP9 complex; (B) Danshenxinkun D-AKT1 complex; (C) 2-isopropyl-8-methylphenanthrene-3,4-dione-AKT1 complex; (D) Sesamin-AKT1 complex; (E) Luteolin-MMP9 complex. In each panel, the compound is shown in blue, the protein in purple, and the interacting protein residues in pink. Please click here to view a larger version of this figure.

Figure 3: Molecular Dynamics (MD) simulations of the main targets and active components of QUF3. (A) Changes in the root-mean-square deviation (RMSD) of the protein over time. (B) Changes in the solvent-accessible surface area (SASA) of protein over time. (C) Changes in the radius of gyration (Rg) of protein over time. (D) Changes in the hydrogen bonds of the protein over time. Please click here to view a larger version of this figure.

Figure 4: Principal Component Analysis (PCA), Free Energy Landscape (FEL), and Dynamic Cross-Correlation Matrix (DCCM) of the main targets and active components of QUF3. (A) PCA of the main targets and active components of QUF3. (B) FEL of the main targets and active components of QUF3. (C) DCCM of the main targets and active components of QUF3. Please click here to view a larger version of this figure.

Table 1: The vital ingredients of QUF3 in the network diagram. 
List of key active constituents identified from the herb-constituent-target network, including compound names, Degree values, and corresponding target proteins. Please click here to view a larger version of this Table.

Energy (KJ/mol)	MMP9-Sesamin	MMP9-Luteolin	AKT1-Sesamin	AKT1-Danshenxinkun D	AKT1-2-isopropyl-8-methylphenanthrene-3,4-dione
Van der Waals Energy	−243.764	−153.616	−211.573	−172.008	−117.13
Electrostatic energy	−4.665	−13.67	−14.797	−8.53	−13.413
Polar solvation energy	82.397	135.935	158.627	122.437	98.14
Nonpolar solvation Energy	−18.82	−18.949	−23.292	−20.987	−16.453
Total Binding Energy	−184.853	−50.3	−91.035	−79.089	−48.856

Table 2: MMPBSA analysis of protein and ligand.

Binding free energies (in kJ/mol) and their components (van der Waals, electrostatic, polar solvation, and non-polar solvation) for the five top-ranked protein-ligand complexes.

Supplementary Table 1: The active ingredients of QUF3.
List of active compounds identified from each herb in QUF3, including Mol ID, compound name, oral bioavailability (OB, %), drug-likeness (DL) value, and the herb of origin. Compounds were screened with OB ≥30% and DL ≥ 0.18. Please click here to download this file.

Supplementary Table 2: The information of ten core proteins
PDB codes, predicted pocket coordinates (X, Y, Z in Å), and grid box dimensions (in Å) were used for molecular docking of the ten core target proteins. Please click here to download this file.

Discussion

There is growing evidence that the endometrium of patients with PCOS exhibits impaired function, which may correlate with higher rates of implantation failure and adverse pregnancy outcomes⁴⁰. Endometrial differences in women with PCOS are mainly characterized by decreased expression of pinopodes, nucleolar channel systems, estrogen receptors, and progesterone receptors during the window of implantation (WOI)⁴¹. Moreover, the expressions of homeobox A10 (HOXA-10), homeobox A11 (HOXA-11), integrin alpha V beta 3 (αvβ3), and leukemia inhibitory factor (LIF) endometrial markers are also reduced, which may lead to infertility in PCOS women^42,43,44. Myo-inositol (MYO) is a novel uterine endometrium marker related to PCOS. As an insulin sensitizer, MYO participates in insulin signal transduction and activates adenosine monophosphate-activated protein kinase (AMPK), thereby inhibiting the decrease in AMPK levels in the uterine endometrium of PCOS women^45,46,47.

Endocrine abnormalities such as hyperandrogenemia (HA), insulin resistance (IR), and obesity are often present in women with PCOS. HA can affect genes and proteins (HOXA, aVβ3 integrin, MECA-79, MAGEA-11, and CDK signaling pathway) involved in endometrial receptivity, altering the endometrial environment and thereby making it unfavorable for embryo implantation⁴⁸. And IR in the endometrium can reduce insulin receptor and glucose transporter levels, causing metabolic abnormalities⁴⁹. Moreover, chronic inflammatory factors and fat factors (such as leptin) further exacerbate the dysfunction of the endometrium, thereby increasing the risk of miscarriage^50,51.

The preferred treatment for endometrial receptivity of PCOS includes lifestyle modifications (dietary management and exercise), drug therapy (metformin), ART, and so on. However, there are still limitations. A large number of studies have confirmed that TCM can effectively improve endometrial receptivity, thereby further improving pregnancy outcomes in PCOS^52,53,54. The occurrence of reduced endometrial receptivity in PCOS is mainly due to kidney deficiency and the obstruction of phlegm and dampness. The QUF3 consists of 7 traditional Chinese herbs that work together to tonify the kidneys, enrich the essence, nourish the blood, and promote blood circulation. In clinical practice, QUF3 has been shown to improve endometrial receptivity in patients with PCOS, but the specific mechanism of action remains unclear. Therefore, using the network pharmacology approach, this study delved into the molecular mechanism by which QUF3 improves endometrial receptivity in PCOS at the microscopic level.

Unlike the traditional “single component–single target” reductionist approach commonly used in Chinese herbal medicine research, network pharmacology can analyze, at a systemic level, the holistic “multi-component, multi-target, multi-pathway” regulatory characteristics of formulations such as QUF3. It helps identify key active components and core targets, thereby narrowing the scope for subsequent experimental validation. However, this method relies on the completeness of annotations in existing databases and cannot simulate the dynamic in vivo metabolic processes of drugs. Static molecular docking can evaluate binding affinity but fails to reflect the conformational changes of protein-ligand complexes over physiological timescales. Molecular dynamics (MD) simulations address this limitation: by calculating RMSD, Rg, SASA, hydrogen bond dynamics, and binding free energy, they enable direct assessment of the stability and flexibility of target protein–ligand complexes, providing trajectory‑based dynamic evidence for judging “strong binding.” Combining network pharmacology with MD simulations retains the efficiency of high-throughput screening while adding dynamic validation of key binding events, forming a cost-effective “broad-first, refined-second, computation-driven” paradigm for investigating mechanisms of TCM formulas. This workflow is highly extensible: replacing the active ingredient library and disease target list allows it to be adapted to other TCM formulas or natural products. All software used is open‑source or freely available, and the key parameters and checkpoints have been clearly described in the Protocol section, facilitating rapid replication by other researchers⁵⁵. The key steps of this protocol include network pharmacology screening and molecular dynamics simulation. Network pharmacology sets a molecular docking scoring threshold of ≤ −5.0 kcal/mol; molecular dynamics simulation uses convergence criteria that require the RMSD trajectory to reach a plateau with stable fluctuations and for key hydrogen bonds to maintain high occupancy throughout the simulation. The key checkpoints are that the docking binding energy meets the threshold, the RMSD trajectory is stable, and the binding free energy calculation passes validation. For those with limited computational resources, the MD simulation time can be appropriately shortened after confirming RMSD stabilization, though this may affect sampling of rare conformations.

As natural flavonoid compounds, quercetin and kaempferol can regulate glucose and lipid metabolism through multiple pathways^56,57. Simultaneously, quercetin and luteolin enhance the endometrial decidualization in mice through anti-aging, anti-inflammatory, and antioxidant stress properties, and improve the pregnancy outcome^58,59. Sesamin has attracted increasing attention due to its immunomodulatory and anti-inflammatory properties⁶⁰. Moreover, sesamin also exhibits the activity of plant estrogen, which can significantly increase the expression of estrogen receptor β in the uterus of rats⁶¹. And β-sitosterol regulates endometrial receptivity in PCOS by modulating the intestinal microbiota and balancing hormonal levels. As the key constituents of Radix Salviae (Danshen), tanshinone IIA, Cryptotanshinone, and danshenxinkun d have anti-inflammatory properties⁶², and tanshinone IIA enhances endometrial receptivity and resolves embryo implantation disorders by inhibiting the glucocorticoids signaling pathway and promoting angiogenesis⁶³.

The PPI network showed that the active components in QUF3 may function through the core targets such as AKT1, EGFR, TNF, TP53, IL6, BCL2, ESR1, IL1B, STAT3, and MMP9. AKT1 is a downstream effector molecule of the PI3K-Akt signaling pathway. Its activated state can trigger a series of cellular biological responses, mainly including cell survival, proliferation, and metabolic regulation^64,65,66. Furthermore, AKT1 also maintains the integrity of the endometrium by reducing oxidative stress and inflammatory responses, which play a crucial role in the development of the endometrium and the maintenance of fertility⁶⁷. EGFR is an important signaling molecule in cell signaling; it is involved in biological processes such as endometrial proliferation and decidualization⁶⁸. The subtype of TNF, TNF-a is closely related to endometrial receptivity in patients with PCOS-induced infertility⁶⁹. An increasing number of studies have shown that TP53 is crucial for embryonic implantation through transcriptional up-regulation of uterine LIF⁷⁰. IL-6 is a pleiotropic cytokine and is responsible for many physiological processes⁷¹. Within the uterine environment, IL-6 triggers the janus kinase JAK-STAT pathway, resulting in STAT3 phosphorylation, whose activation is necessary for implantation^72,73. BCL2 is an important anti-apoptotic gene that maintains the survival of endometrial cells by inhibiting the apoptotic pathway⁷⁴. ESR1 can activate downstream signaling pathways such as the PI3K-Akt signaling pathway, promoting the survival, proliferation, and differentiation of endometrial cells⁷⁵. Studies have shown that MMP9 is the main trigger of blastocyst implantation. MMP9 is typically expressed in the epithelial cells of the endometrium and plays a significant role in endometrial receptivity⁷⁶.

The KEGG pathway enrichment analysis of QUF3 revealed that the pathways mainly involved multiple signaling axes including PI3K-Akt, MAPK, HIF (HIF-1 and HIF-2), FoxO, JAK-STAT, and Ras/Rap1, as well as cellular senescence, focal adhesion, thyroid hormone and relaxin signaling, cytokine-cytokine receptor interaction, and immune-related pathways such as IL-17, TNF, T cell receptor, and C-type lectin receptor signaling, along with osteoclast differentiation and prolactin signaling, which may involve anti-inflammatory effects, antioxidative stress, regulation of immune responses, inhibition of cell apoptosis, and enhancement of endometrial response⁷⁷. MD revealed robust binding activity between the core targets and their respective active components, thereby corroborating the network pharmacology predictions. Subsequent MD simulations, assessed through RMSD, SASA, Rg, principal component analysis, and free energy landscape analyses, confirmed that the five top-ranking complexes from the docking studies exhibited sufficient dynamic stability and flexibility. Their protein-ligand interactions were primarily governed by hydrogen bonds and van der Waals forces.

This study still has several limitations. First, network pharmacology relies on existing database information, which may not fully account for in vivo metabolism or complex regulatory networks, potentially introducing biases. Second, MD simulations were performed only for selected single targets and may not fully capture the multi-component synergistic effects of TCM. Third, the binding energy threshold (−5.0 kcal/mol) for “strong activity” is empirical; different thresholds may yield different rankings. Finally, all findings are computational predictions and require further in vitro, in vivo, and clinical validation. Future work should include experimental verification of the predicted mechanisms and expand MD simulations to additional target-ligand pairs to better reflect the holistic effects of QUF3.

The mechanism of action of QUF3 in treating endometrial receptivity in PCOS may be associated with components such as luteolin, sesamin, and danshenxinkun d in the formula, which may provide a reference for further experimental verification.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (Version 2.3)	Lab of Systems Pharmacology	https://www.tcmsp-e.com/load_intro.php?id=43	Database for screening of active components
Swiss Target Prediction online database (2019 version)	Molecular Modelling Group, University of Lausanne & SIB Swiss Institute of Bioinformatics	http://www.swisstargetprediction.ch/	Online target prediction database
PubChem database	National Library of Medicine, NIH	https://pubchem.ncbi.nlm.nih.gov/	Chemical database for compound SMILES/SDF retrieval
UniProt database	NIH	https://www.uniprot.org/uniprotkb	Protein database for UniProt IDs and gene names
Online Mendelian Inheritance in Man (OMIM)	Johns Hopkins University	https://www.omim.org/	Disease target database
DrugBank	University of Alberta	https://www.drugbank.ca/	Disease target database
Therapeutic Target Database (TTD)	Zhejiang University	https://ttd.idrblab.cn/	Disease target database
GeneCards (Version 5.26)	Weizmann Institute of Science, LifeMap Sciences	https://www.genecards.org/	Disease target database (relevance score >5)
Venny visualization platform (Version 2.1.0)	Centro Nacional de Biotecnología (CNB-CSIC)	https://bioinfogp.cnb.csic.es/tools/venny/	Online Venn diagram tool for intersection target visualization
Cytoscape (Version 3.10.3)	National Resource for Network Biology, NHGRI	https://cytoscape.org/	Network analysis and visualization software for H-C-T, PPI, and pathway networks
STRING database (Version 12.0)	STRING Consortium	https://string-db.org/	Protein-protein interaction database for PPI network construction
DAVID database (Version v2025-2)	U.S. Department of Health & Human Services, NIH	https://davidbioinformatics.nih.gov/	Gene enrichment analysis database for GO and KEGG
AutoDockTools (Version 1.5.7)	The Scripps Research Institute	https://autodock.scripps.edu/	Molecular docking software for preprocessing and calculation
PyMOL software (Version 3.0.3)	Schrödinger	https://pymol.org/	Molecular visualization software for docking visualization
Protein Data Bank (PDB)	RCSB	https://www.rcsb.org/	Protein structure database for downloading target PDB files
GROMACS (Version 2023.2)	GROMACS development team	https://www.gromacs.org/	Molecular dynamics simulation software
VMD (Version 1.9.3)	University of Illinois	https://www.ks.uiuc.edu/Research/vmd/	Molecular visualization and trajectory analysis software
gmx_MMPBSA	Open-source (GitHub)	https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA	Binding free energy calculation tool from MD trajectories
Science Data Bank	Chinese Academy of Sciences, Computer Network Information Center	https://www.scidb.cn/	Scientific data repository for data availability and deposition