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 Officinal...

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 ...

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