Piperine’s potential in treating polycystic ovarian syndrome explored through in-silico docking

Polycystic Ovarian Syndrome (PCOS) is a multifaceted metabolic and hormonal condition that impacts women in their procreative ages, identified by ovarian dysfunction, hyperandrogenaemia overweight and insulin insensitivity. The piperine, an important alkaloid compound of black pepper has shown promise in modulating various physiological processes. In this work, employed computational docking studies to explore the potential of piperine as a treatment for PCOS. Utilizing computational methods, we analyzed the binding interactions between piperine and key molecular targets implicated in PCOS pathogenesis, including hyperandrogenism, and “oligomenorrhea. The network pharmacology analysis report found 988 PCOS-related genes, 108 hyperandrogenism-related genes, and 377 oligomenorrhea-related genes, and we finally shortlisted 5 common genes in PCOS, hyperandrogenism, and “oligomenorrhea”: NR3C1, PPARG, FOS, CYP17A1, and H6PD. Our results reveal favorable binding affinities with PPARG (-8.34 Kcal/mol) and H6PD (-8.70 Kcal/mol) and interaction patterns, suggesting the potential of piperine to modulate these targets. Moreover, the reliability of the piperine-target interactions was revealed by molecular simulations studies. These findings support further experimental investigations to validate the therapeutic efficacy of piperine in PCOS management. The integration of computational approaches with experimental studies has the potential to lay the groundwork for the creation of new therapies specifically targeting PCOS and related endocrine disorders. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-024-72800-6.


Identification of key drug targets in PCOS by network pharmacology approach
The DisGeNET database(https://www.disgenet.org/home/)was queried for genes connected with polycystic ovarian syndrome (PCOS), using "Polycystic Ovarian disease, " "Hyperandrogenomia, " and "Oligoanovalation" as search keywords.The outcomes of the three data sets were combined and de-duplicated to get genes that particularly matched the three search criteria.

Construction of disease target genes analysis
The study employed the online Venn diagram application at (https://bioinfogp.cnb.csic.es/tools/venny/)(Oliveros, J.C. (2007-2015)to illustrate the points where disease hits meet 13 .These potential therapeutic targets were subsequently utilized to construct the KEGG Path Examination, Gene Ontology Predicting, and Protein-Protein Interacting Networks.

Constructing interactions network of protein-protein
Using STRING (https://cn.string-db.org/) the intersection gene was imported and establishes the interaction connection between target proteins in order to identify the network's primary targets 14 .Constructing a network of protein-protein interactions was built by cytoscape 15 .The degree value was determined and the core targets were selected via Cytoscape's own analysis.

Molecular docking study
Approach for carrying out molecular docking investigations with an adjusted flexible docking methodology by Rizvi et al. (2013)  21 .It involves using the MGL graphic tool with AutoDock Tool 4.2 Version 22 for exploring interactions of piperine with target proteins shown in Table 1.Proteins are prepared from PDB files, and receptor grids are created, adapting to binding pockets.Docking parameters include allowing ligand rotation and selecting optimal docking postures based on RMSD, Ki, and binding energies.Cygwin 3.5.3.Version 23 software is used for manual comparison, with ten configurations generated per protein-ligand combination.An exhaustiveness of 10 is used, and Discovery Studio 2017 is employed for post-docking analysis of ligand-protein interactions 21 .

Molecular simulation study
The molecular simulation study applied to the docking study-identified ligand with the greatest binding ability, using AutoDock software 24 .Hits were analyzed using the Desmond 2.2 Version Program for Molecular Dynamics 25 .Each protein-ligand combination was enclosed in an orthorhombic box and hydrated using TIP3P water molecules.On addition, 0.15 M Sodium ions (Na+) and chloride ions (Cl-) were introduced into the system to neutralize it using the OPLS3e force field.The standard Desmond equilibration protocol was followed, including steps such as Brownian Dynamics NVT simulation, NPT simulation with constraints, and a final NPT simulation without constraints, each at a temperature of 10 K.The final MD simulation ran for 100 nanoseconds at standard temperature (300 K) and pressure (1.013 bars) 25 .Pressure and temperature both controlled utilizing the Nose-Hoover Chain thermometer and Martyna-Tobias-Klein barometer 26,27 .Using a smooth Particle Mesh Ewald (PME) approach with a RESPA integrator and a time step of two femtoseconds, long-distance electrostatic forces were calculated.Every 10 ps, coordinates and energy were saved for trajectory analysis.

Results
The 3D structure of piperine was depicted in Fig. S1a.Table 2 presents the physiochemical properties of piperine, while Table S1 summarizes its ADMET properties.The ADMET properties shown that piperine has no toxicity at end points in silico models.The pkcsm tools provided predictions for a range of properties concerning piperine, including adsorption, distribution, release, waste elimination, and toxic effects represents in Table S1 28,29 .

Gene ontology's prediction
In fold enhancement analysis, the False Detection Ratio (FDR) is calculated using the nominal P-value obtained from the hyper geometric test.In order to calculate fold enrichment, divide the proportion of your list's genes that belong to a certain pathway by the proportion of background genes in that route.FDR reveals the probability of the enrichment occurring by chance.Large paths often have reduced FDRs because of improved statistical  power.Fold Enrichment chart (Fig. 1a) shows how significantly overrepresented a certain pathway's genes are as a metric of impact magnitude.The charts for biological activities (Fig. 1b), molecular mechanisms (Fig. 1c), and cell component (Fig. 1d) showed the projected GO keywords.

Protein-protein interaction network analysis
We utilized PPI networks to investigate the connections that exist between various gene targets and to locate important network genes.With a confidence level of > 0.900, Homo sapiens was the chosen species, and five common targets of proteins were inserted into STRING, (Fig. S1c) shows that there were 35 nodes, 288 edges, a 16.5 typical node degree, a 0.781 the typical local clustered ratio, 87 expected edges, and a p-value of < 1.0e-16 for PPI enrichment 27 .Protein interactions involving Nuclear Receptor Subfamily 3 Group C Member 1 (NR3C1), Peroxisome Proliferator Activated Receptor Gamma (PPARG), Transcription Factor AP-1 Subunit C-Fos (FOS), Cytochrome P450 Family 17 Subfamily A Member 1 (CYP17A1), and Hexose-6-Phosphate Dehydrogenase/ Glucose 1-Dehydrogenase (H6PD) were directly observed 30 .So, we predicted 5 common genes in "Poly Cystic   Ovary Syndrome (PCOS)", "Hyperandrogenism" and "Oligomenorrhea".The five genes that are cross-targets to the piperine for the treatment of PCOS.

KEGG enrichment pathway analysis
Several KEGG pathways were found to be closely linked with PCOS, including Prolactin signaling, Pentose phosphate pathway, Osteoclast differentiation, Thyroid cancer, endocrinal hormones biosynthesis, Non-alcohol fatty liver disease, Ovarian steroidogenesis, Lipid metabolisms and atherosclerosis, Cortisol synthesis and secretion, and Amphetamine addiction.Table S2 illustrates the 30 pathways deemed most crucial for PCOS 31 .Furthermore, Fig. 1e displays the top pathways with statistically significant differences (P < 0.05) as identified by the KEGG enrichment analysis.

Molecular docking analysis
Following extensive validation utilizing the Lipinski rule of five targets and ADMET properties, we performed docking of piperine against the target proteins, including glucose 1-dehydrogenase/hexose-6-phosphate dehydrogenase (H6PD), Subunit C-Fos of  3. The docked complexes of all the protein targets were shown in (Fig. 2a-e) 32 .The 2D interaction analysis of targeted proteins with the piperine we analyzed the  protein-ligand docked complexes shown in (Fig. 3a-e).This suggests that target binding necessitates the amino acid residues.

Molecular dynamics (MD) simulations
The Protein Data Bank (PDB) typically houses structural data of the target biomolecule obtained through X-ray or NMR techniques.Computational modeling approaches, such as homology strategies, are commonly employed to derive coordinate and geometric information regarding protein structure 33 .This step also considers the molecular environment, including solvation and ionic strength.The Maxwell-Boltzmann distributions at the optimal simulation temperature govern the initial velocities of molecules 34 .Interactions between atoms are computed using a force field, which defines both non-bonded and bonded terms based on the atomic composition, including electrostatic, Lennard-Jones, dihedral, angle, improper dihedral, and bond terms 35 .
The overall potential energy of the system is determined by accounting for all potential interactions between atoms within the protein system.
The force exerted is determined by deriving the potential energy function, equating to acceleration divided by mass according to Newton's second law.This derivative yields a set of 3 N coupled second-order ordinary differential equations that must be solved mathematically, despite inherent complexities.The solution involves advancing positions and velocities using a numerical algorithm, iterating through time steps to generate molecular dynamics (MD) trajectories with constant energy.Coupling the system to a thermostat and employing established methods like Langevin dynamics or Nose-Hoover algorithms ensures consistent temperature dynamics 35 .Several prominent force fields, such as CHARMM (www.charmm.org),AMBER (www.ambermd.org),and GROMOS (www.gromacs.org),along with molecular dynamics simulation packages like NAMD (www.ks.uiuc.edu/Exploration/namd/)and VMD (www.ks.uiuc.edu/Exploration/vmd/),significantly enhance macromolecular MD simulations 36 .Analyzing statistical properties under various initial and external conditions, such as hydrogen bond networks, hydrophobic interactions, and solvent accessible surface area, is achievable through examination of molecular dynamics trajectories 37 .By assessing the spatial relationships between hydrogen bond donors and acceptors with fixed cutoff angles and bond lengths, it becomes feasible to identify stable hydrogen bonds from transient interactions 38 .Additionally, quantifying hydrophobic stabilization effects can be achieved through solvent-accessible surface area analysis, which involves delineating surface regions using a 1.4 Å probe radius and calculating statistical metrics 39 .
Empirical energy functions, employed in creating force fields, are customizable and tailored to specific classes of molecules through parameterization.However, their applicability is limited to native systems and environments, making them unsuitable for modeling non-native systems.As a result, results often necessitate comparison with experimental data to validate accuracy and identify areas for methodological improvement.This ongoing process underscores the continuous advancement of foundational force fields and methodologies 38 .
An additional critical consideration is the ability to adequately sample the vast array of possible conformations available to even the simplest biomolecules.One potential limitation of molecular dynamics simulations is the typically restricted maximum time step usable for simulation (around 2 femtoseconds).Protein monomers' solved structures typically encompass 40,000 atoms, while structures of dimers and membrane-bound proteins can range from 200,000 to 500,000 atoms.Achieving simulation times in the microsecond range and beyond for such large proteins with current hardware and software requires algorithmic improvements and leveraging high-performance computing resources 40 .

Root mean square deviation (RMSD)
The Root Mean Square Deviation (RMSD) of the protein and ligand from their reference positions in the bound complex, achieved through optimal superimposition of receptor structures, is widely accepted for assessing the accuracy of docking geometries 41 .The protein's all-atom RMSD provides insight into how ligand binding affects the overall 3D structure.Initially, all structures were aligned based on non-hydrogen atoms, and RMSD values indicate changes in protein structure relative to a reference state.Comparisons were made with initial and average structures obtained from simulations 41 .An RMSD value ≤ 2 Å is considered optimal and recommended 42 .Interestingly, the lowest energy conformer for each ligand did not always exhibit the best pose (i.e., lowest RMSD).In some cases, a second conformer with a slightly higher binding affinity (0.0-0.5 kcal/mol) displayed a lower RMSD than the lowest energy conformer.
Initially, during a 100 ns molecular dynamics simulation, trajectories of a protein or receptor, whether in its apo form or complexed with ligands, were analyzed to assess the stability of these proteins, ligands, and their molecular interactions.A comparison was conducted between the MD trajectories of the apoprotein and the ligand-bound complex to evaluate protein fluctuations (RMSD).The graphical representation of RMSD values obtained for each ligand's predicted lowest binding energy models were illustrated in Fig. S2.

Protein RMSD
The plot illustrates the RMSD evolution of a protein (left Y-axis).All protein frames are initially aligned to the reference frame backbone, and RMSD is subsequently calculated based on selected atoms.The RMSD Monitoring protein provides insights into its structural conformation throughout the simulation.RMSD analysis helps determine if the simulation has equilibrated fluctuations toward the end of the simulation should reflect a thermal average structure.Changes on the order of 1-3 Å are generally acceptable for small, globular proteins.Larger changes indicate significant conformational shifts during the simulation.
Additionally, it's crucial for the simulation to converge -RMSD values should stabilize around a fixed value.If protein RMSD continues to increase or decrease towards the end of the simulation, it suggests that equilibration hasn't been achieved, and the simulation may require additional time for rigorous analysis.

Ligand RMSD
The plot above displays the RMSD of a ligand (right Y-axis), indicating its stability relative to the protein and its binding pocket.In 'Lig fit Prot' , the ligand's RMSD is computed after aligning the protein-ligand complex to the reference protein backbone, followed by measuring the RMSD of the ligand's heavy atoms.If the observed values are notably larger than the RMSD of the protein, it suggests that the ligand may have moved away from its initial binding site.This analysis helps assess the stability of the ligand within the protein's environment during the simulation.
The RMSD plot for each of the four ligands is depicted in Fig. S2.In this study, RMSD values were calculated specifically for the C-alpha atoms of the protein structures.The RMSD values of the four ligands and those of the protein were compared within each complex.An acceptable and stable deviation between the RMSD of the protein alone and the RMSD of the protein-ligand complex typically ranges from one to two angstroms.RMSD serves as a measure of protein stability throughout the simulation period.

Root mean square fluctuation (RMSF)
Following roto-translational least-squares fitting, the global structural similarity of commonly used macromolecules (proteins) is assessed by their root mean square deviation (RMSD).Conversely, experimental x-ray B-factors are utilized to focus on structural heterogeneity at the atomic level in macromolecules, providing direct information about root mean square fluctuations (RMSF) that can also be computed from molecular dynamics simulations.To evaluate the flexibility and dynamics of amino acids over a 100 ns simulation, Root Mean Square Fluctuation (RMSF) values were calculated.Each amino acid residue's RMSF value is examined to understand its interaction dynamics with the ligand 43 .
The C-alpha backbone of the protein was analyzed based on RMSF values for each amino acid residue to further investigate the flexibility of individual residues and their interaction dynamics in ligand-bound and ligand-unbound simulations.A computational RMSF prediction indicating negative values suggests reduced fluctuations in the presence of ligands or increased fluctuations in their absence 44 .This analysis helps characterize how ligand binding influences the structural dynamics of the protein at the atomic level.

Protein RMSF calculation
The root mean square fluctuation (RMSF) calculation provides insight into the residual mobility of amino acids.Amino acids with higher RMSF values typically indicate greater flexibility and fluctuation within the receptor structure.Conversely, those with lower RMSF values are indicative of more rigid and stable regions within the receptor.
The RMSF values of the docked complex (Hexose-6-phosphate dehydrogenase-Piperine) have been calculated and are depicted in Fig. S3a.The Fig. S3 illustrates how different amino acids within the complex exhibit varying levels of mobility, offering a detailed view of the receptor's dynamic behavior in response to ligand binding.
When comparing to the native protein structure, the average RMSF values of receptor amino acid residues tend to be lower when bound to bioactive compounds like piperine and Hexose-6-phosphate dehydrogenase.This observation suggests that the docked ligand has contributed to the formation of a rigid and stable complex.This stability in RMSF values indicates that the ligand binding has ind uced a more constrained and structured conformation in the receptor, potentially influencing its functional behavior.
The RMSF plot of Hexose-6-phosphate dehydrogenase docked with the ligand piperine is depicted in Fig. 4.During the 100 ns MD simulation study, amino acids in the protein structure exhibited varying RMSF values, indicating fluctuation and rigidity.Specifically, residues in the ranges of 60 to 75, 250 to 300, 350 to 380, and 410 to 430 showed the highest fluctuations during the simulation.These residues are predominantly located on the cytoplasmic side of the C-terminal helix.In contrast, other regions of the protein displayed lower fluctuations, with certain areas, such as the C-terminal helix, showing less variability as depicted in the plot.
The Protein Root Mean Square Fluctuation (L-RMSF) is useful for characterizing changes in the protein and ligand atom positions.The RMSF for atom i is: On the RMSF plot, peaks indicate areas of the protein that exhibit the highest fluctuations during the 100 ns simulation study.Typically, the N-and C-terminal tails fluctuate more than other regions of the protein.
Secondary structure elements such as alpha helices and beta strands are generally more rigid and thus exhibit lower fluctuations compared to the unstructured loop regions, which tend to be more flexible and dynamic.

Ligand RMSF calculation
The atom-by-atom breakdown of the ligand's fluctuations in the Ligand RMSF plot corresponds to the 2D structure depicted in the top panel.This analysis provides insights into the entropic contributions of ligand fragments during the binding event and their interactions with the protein.
In the bottom panel, the 'Fit Ligand on Protein' line shows the ligand fluctuations with respect to the protein.The protein-ligand complex is first aligned based on the protein backbone, and then the ligand RMSF is measured for the ligand's heavy atoms (Fig. S3b).This alignment allows for a detailed examination of how each part of the ligand behaves in the context of the bound complex.
The Ligand Root Mean Square Fluctuation (L-RMSF) is uses to characterize the changes present in the protein and ligand atom positions.

RM SF
where T is the trajectory time over which the RMSF is calculated, t ref is the reference time, r i is the position of residue i; r' is the position of atoms in residue i after superposition on the reference, and the angle brackets indicate that the average of the square distance is taken over the selection of atoms in the residue.

Protein-ligand contact residues
Protein-ligand interactions can be tracked during the simulation and are classified by type, as illustrated in Fig. 5a.Hydrogen bonds, hydrophobic contacts, ionic interactions, and water bridges are the four basic types of these interactions.The 'Simulation Interactions Diagram' panel can be utilized to investigate the different subtypes that are present in each category in more detail.Throughout the simulation trajectory, the stacked bar charts in the graphic are normalized.When an interaction has a value of 0.7, for example, it means that it continues for 70% of the simulation.Since some protein residues may make several interactions of the same subtype with the ligand, values greater than 1.0 may occur.
This detailed categorization and visualization help in understanding the stability and nature of protein-ligand interactions, providing insights into the binding mechanisms and dynamics of the complex.
Hydrogen bond contacts H-bonds, or hydrogen bonds, are essential for ligand binding.Because of their profound effects on drug selectivity, metabolism, and absorption, hydrogen-bonding qualities must be taken into account while designing new medications.Four subtypes of hydrogen bonds exist between a ligand and a protein: side-chain donor, sidechain acceptor, and backbone acceptor.The present geometric requirements for protein-ligand H-bonds are as follows: a donor angle of ≥ 120° between the donor-hydrogen-acceptor atoms (D-H•••A); an acceptor angle of ≥ 90° between the hydrogen-acceptor-bonded atom atoms (H•••A-X); and a distance of 2.

Å between the donor and acceptor atoms (D-H•••A).
Hydrogen bonding, a type of polar binding, occurs when a hydrogen atom is bound to a highly electronegative molecule and engages in electrostatic interactions with a hydrogen bond acceptor 45 .Hydrogen bonds (H-bonds) are crucial stabilizing interactions between two molecules at the binding site of a protein.To examine these interactions, the number of H-bonds is determined 46 .The H-bonds were recorded throughout the 100 ns of the receptor-ligand MD simulations (Fig. 5a).
A hydrogen bond network is a chain of hydrogen bonds connecting the side chains of multiple residues.The side chain molecules that can form hydrogen bonds are displayed in the Sidechain interaction.Molecular Dynamics (MD) simulations have also been used to explore protein hydration 47 .Hydrogen bonding is one of the most critical interactions in a biological system for stabilizing the protein-ligand complex.Figure 5b depicts the number of hydrogen bonds formed in protein-ligand complexes over the last 100 ns.The hydrogen bonding network demonstrates that these bonds are steady interactors for proteins.
In the case of the ligand (piperine) complexed with Hexose-6-phosphate dehydrogenase (H6PD), the number of intermolecular H-bonds is greater during the last 100 ns of simulation, indicating stable interaction and binding.

Hydrophobic contacts
Hydrophobic contacts fall into three subtypes: π-Cation; π-π; and Other, non-specific interactions.Generally these type of interactions involve a hydrophobic amino acid and an aromatic or aliphatic group on the ligand, but we have extended this category to also include π-Cation interactions.The current geometric criteria for hydrophobic interactions is as follows: π-Cation -Aromatic and charged groups within 4.5 Å; π-π -Two aromatic groups stacked face-to-face or face-to-edge; Other -A non-specific hydrophobic sidechain within 3.6 Å of a ligand's aromatic or aliphatic carbons.

Ionic interaction
Ionic interactions or polar interactions, are between two oppositely charged atoms that are within 3.7 Å of each other and do not involve a hydrogen bond.We also monitor protein-metal ligand interactions, which are defined by a metal ion coordinated within 3.4 Å of protein's and ligand's heavy atoms (except carbon).All ionic interactions are broken down into two subtypes: those mediated by a protein backbone or side chains.

Water bridges
Water Bridges are hydrogen-bonded protein-ligand interactions mediated by a water molecule.The hydrogenbond geometry is slightly relaxed from the standard H-bond definition.The current geometric criteria for a protein-water or water-ligand H-bond are: a distance of 2.8 Å between the donor and acceptor atoms (D-H•••A); a donor angle of ≥ 110° between the donor-hydrogen-acceptor atoms (D-H•••A); and an acceptor angle of ≥ 90° between the hydrogen-acceptor-bonded atom atoms (H•••A-X).

Protein-ligand contacts timeline
A timeline representation of the interactions and contacts (H-bonds, Hydrophobic, Ionic, Water bridges) were summarized in Fig. 5a.The top panel shows the total number of specific contacts the protein makes with the ligand over the course of the trajectory.The bottom panel shows which residues interact with the ligand in each trajectory frame.Some residues make more than one specific contact with the ligand, which is represented by a darker shade of orange, according to the scale to the right of the plot (Fig. 5b).

Ligand torsion profile
Each ligand in this study has only one rotatable bond.In our 100 ns MD simulation, the ligand (Piperine), a synthetic compound, moves without breaking its bonds.One of the movements observed is rotamerization, which involves rotations occurring around single bonds in the ligand 48 .The ligand torsion profile provides data on how well the ligand torsion aligns with theoretical estimations and indicates the preferred positions of the ligand throughout the simulation 49 .The ligand torsion profile plot is represented in Fig. S4.
This plot helps in understanding the flexibility and conformational preferences of the ligand, which are crucial for its binding affinity and specificity.By analyzing the torsion profile, we can gain insights into the stability and dynamics of the ligand within the binding site of the protein.
Plots with dials, also called radials, show the torsion conformation during the simulation.The time progression is shown radially outward from the center of the radial plot, where the simulation begins (Fig. S4).
By showing the probability density of the torsion, the bar plots give an overview of the data from the dial plot.Plotting the potential of the corresponding torsions sums up to represent the potential of the rotatable bond when torsional potential information is provided.The potential values are displayed on the chart's left Y-axis and are given in kcal/mol.The conformational strain that the ligand undergoes in order to retain a protein-bound shape can be understood by looking at the histogram and torsion potential correlations (Fig. S4).

Ligand (piperine) properties
I. Ligand RMSD: Root mean square deviation of a ligand with respect to the reference conformation (typically the first frame is used as the reference and it is regarded as time t = 0).II.Radius of Gyration (RoG): Measures the 'extendedness' of a ligand, and is equivalent to its principal moment of inertia.We have determined the RoG in light of the characteristic elements of apoprotein-ligand buildings (Fig. S5a).Over the 100 ns MD simulation that was investigated, all of the complexes have a constant average RoG.Such outcomes demonstrate that greater part of ligands have framed conservative and stable interaction with apoprotein when contrasted with the reference protein structure.III.Intramolecular Hydrogen Bonds (intraHB): Number of internal hydrogen bonds (HB) within a ligand molecule.IV.Molecular Surface Area (MolSA): Molecular surface calculation with 1.4 Å probe radius.This value is equivalent to a van der Waals surface area.V. Solvent Accessible Surface Area (SASA): Surface area of a molecule accessible by a water molecule.Solvent accessible surface area (SASA) of proteins has everlastingly been considered as a conclusive Fig. S5a protein unfolding and stability evaluation.The surface around a protein is portrayed by the particle's van der Waals contact surface and a speculative focus of a dissolvable sphere 50 .VI. Polar Surface Area (PSA): Solvent accessible surface area in a molecule contributed only by oxygen and nitrogen atoms.

MM/PB(GB)SA calculation
MM/PBSA is a famous strategy to work out calculate the binding affinities in docked complex.In the beginning, it was presented as a method with six well-defined terms and no adjustable parameters.However, since then, a lot of different options have been suggested, and the user now has to make a lot of decisions about things like the dielectric constant, the parameters for the non-polar energy, the PB or GB calculation radii, including or not including the entropy term, MD simulations, or minimizations.In practice, it typically produces results of intermediate quality, typically superior to docking and scoring but inferior to FEP; for instance, r 2 = 0.3 for the entire PDB bind database but r 2 = 0.0-0.8 for individual proteins.The positioning of ligands is frequently unacceptable assuming their affinities contrast by lesser than 12 kJ/mol 51 .Although it allows for greater variation in the ligand scaffold than AP methods and is less sensitive to changes in the net charge, the tested receptorligand system has a significant impact on the results 52 .In our study, the MM/PBSA calculation shows that the docked complex of Hexose-6-phosphate dehydrogenase/Piperine has the binding affinity of -19.71 kcal/mol (Fig. S5b).

Discussion
PCOS is a particularly typical endocrinological illness, identified by persistent anovulation and hyperandrogenism 53 .
If PCOS is not efficiently and adequately treated, it can drastically decrease life expectancy 54 .Many biomarkers have been linked to PCOS and could serve as potential treatment targets.However, the exact mechanism of gene regulation leading to the progression of the syndrome remains unclear 55 .The docking process assesses the binding affinity between molecules.In this research endeavor, we sought to investigate the therapeutic potential of piperine in treating PCOS by employing bioinformatics analysis to gain insights into their biological roles.Piperine, a natural polyphenol derived from black pepper (Piper nigrum) and long peppers (Piper longum), It's been used as a food spice and a potential therapy for several illnesses, including inflammation, obesity, and various malignancies 56 .
Piperine is an orange needle-shaped crystal with melting and lighting points of 131-132 °C.Chloroform and methanol soluble; unable to dissolve in water 57 .Priyanka et al. (2020) found that the results of the ADMET study on piperine derivatives derived from natural sources indicated that most ADME (Absorption, Distribution, Metabolism, and Excretion) features exhibited favorable characteristics, suggesting that these compounds are promising candidates for further development 58 .In the study investigation, the NCBI PubChem database makes the 3D structure of piperine accessible in the format of a structure-data file (SDF).The 3D structure of piperine was seen in Fig. S1a.Using the SwissADME web server, the physiochemical characteristics of piperine were examined, and ADMET characteristics were predicted using the online PKCSM program.Table 2 displayed the physiochemical characteristics of piperine, while Table S1 provided a summary of its ADMET characteristics.In this studiy, we confirmed as same results through in silico ADMET models that piperine exhibits no harmful effects at end points.Kamboj et al. utilized molecular docking to assess the binding affinity between compounds and CYP17.Their findings indicate that all phytocompounds exhibited blocking or hindrance of the CYP17 enzyme's activity, with docked scores ranging from − 3.7 to -9.5 59 .Using computational molecular docking studies, Amudha studied the effectiveness of phytochemicals derived from Cadaba fruticosa (L.) Druce in inhibiting CYP17.Docking scores ranged from − 3.3 to -7.9, indicating that all 20 drugs demonstrated moderate to strong inhibition.In particular, docking scores of -7.9 and − 6.8 indicated strong inhibition of androstan-3-one and 17-hydroxy-2, 4-dimethyl respectively 60 .A more extensive negative value in the docking score denotes a better affinity among the ligand and the protien.The docking of ligands and proteins was ranked according to their binding affinity, with lower (more negative) scores indicating stronger interactions.Upon conducting a computational docking study, we found that piperine has the greatest affinity for binding to the targeted proteins.The docked complexes of all the protein targets are illustrated in Fig. S2 through 13.Over all, the docking study revealed that piperine exhibited the highest binding affinity with Peroxisome Proliferator Activated Receptor Gamma (-8.34 kcal/ mol) 31 and Hexose-6-Phosphate Dehydrogenase/Glucose 1-Dehydrogenase (-8.70 kcal/mol).
The expression of PPAR-γ in the ovaries has been linked to the formation of follicles and ovarian function.Changes in PPAR-γ signalling might interfere with ovarian steroidogenesis and folliculogenesis, which could lead to reproductive problems in PCOS, including irregular menstruation and infertility (Komar, 2005)  61 .In addition, Ahmadian et al. (2013) highlighted the significant function of the PPAR-γ protein in regulating adipocyte development, glucose levels, and lipid homeostasis 62 .An important nuclear receptor recognized as peroxisome proliferator-activated receptor-γ (PPAR-γ) is connected with hyperandrogenemia and plays a role to regulating energy balance, as noted by Stump et al. (2015)  63 .Chen, et al.., (2015) noted that Androgen metabolism and androgen receptor signalling have been associated with PPAR-γ.The hyperandrogenic phenotype associated with PCOS may be attributed to dysregulated PPAR-γ activation 64 .Unluturk et al. (2007)  found that PPAR-γ genes are associated with PCOS occurrences across different ethnic populations.Their study provided evidence linking PPAR-γ with the pathogenesis of PCOS 65 .According to our docking study, piperine exhibits the highest binding affinity with the Gamma-Activated Peroxisome Proliferator Receptor, with a calculated energy of -8.34 kcal/mol.Overall, PPAR-γ docking studies in PCOS research give important insights into the underlying molecular processes of the illness and this could potentially lead to the development of novel treatment approaches.
Docking studies involving Hexose-6-phosphate dehydrogenase (H6PD) genes and their relevance to PCOS are not commonly reported in the literature.While H6PD is involved in various metabolic pathways and play the important role in PCOS pathophysiology, specific docking studies focusing on H6PD genes in the context of PCOS are limited 66 .However, computational docking studies involving other genes or proteins implicated in PCOS, such as PPAR-γ or insulin receptors, have been conducted to explore potential therapeutic interventions.Continuing research exploring the interaction between H6PD genes and pertinent molecular targets in PCOS could yield valuable insights into the fundamental mechanisms of the disorder.
Hexose-6-phosphate dehydrogenase (H6PD) genes play the important role in polycystic ovarian syndrome (PCOS) by influencing various metabolic pathways and hormonal regulation.These genes contribute to glucose metabolism, oxidative stress management, steroid hormone production, insulin sensitivity, and lipid metabolism, all of which are dysregulated in PCOS (Li Yan et al., 2019)  67 .Understanding the involvement of H6PD in PCOS can provide insights into the disorder's mechanisms and potential therapeutic strategies.Cortisone reductase deficiency, caused by inactivating mutations in the enzyme Glucose 1-dehydrogenase/hexose-6-phosphate dehydrogenase (H6PD), mimics of PCOS and manifests as hyperandrogenism that cannot be explained by standard tests (Qin and Rosenfield 2011) 68 .Based on docking study, it was found that piperine demonstrates the greatest binding affinity with the Glucose 1-dehydrogenase/hexose-6-phosphate dehydrogenase (H6PD), showing a calculated energy of -8.70 kcal/mol.
Molecule Dynamic Simulation is regularly utilized to refine proteins or any molecular structures and find previously elusive information about the time evolution of conformations (Vanommeslaeghe et al., 2014)  35 .The protein-ligand RMSD showed that hexose-6-phosphate dehydrogenase docking complexes with 3.6 were present.The RMSD value of piperine when it binds to hexose-6-phosphate dehydrogenase (3.6) is larger than in the natural crystal structure of Glucose 1-dehydrogenase/hexose-6-phosphate dehydrogenase (H6PD).To confirm the adaptability and flexibility of amino acids during simulation for 100 nano seconds, the root of mean squares change (RMSF) was determined.Every residue of the amino acid RMSF value is checked as it cooperates with the ligand along a direction (De Vita et al., 2021) 69 .A well-known method for calculating the binding affinities of docked complexes is MM/PBSA.According to our study's MM/PBSA calculations, the docked complex of hexose-6-phosphate dehydrogenase and piperine exhibits a binding affinity of -19.71 kcal/ mol, as seen in Fig. S3b.

Conclusion
In conclusion, the combination of network pharmacology and molecular docking analyses offers empirical evidence supporting the potential pharmaceutical application of exploring the potential role of piperine in the treatment or management of polycystic ovarian syndrome (PCOS).Through computational analyses, piperine exhibited promising drug-like characteristics and demonstrated favorable interactions with key proteins associated with PCOS pathogenesis, particularly Hexose-6-phosphate dehydrogenase.Molecular dynamics simulations further confirmed the stability and efficacy of the piperine-Hexose-6-phosphate dehydrogenase complex, suggesting its potential as a therapeutic agent for managing PCOS symptoms.These findings lay a solid scientific foundation for further exploration and development of piperine-based treatments for PCOS.

Fig. 1 .
Fig. 1.(a) Bubble chart of Fold enrichment chart of common protein targets in PCOS (b) Biological Process of common protein targets (c) Molecular Function Process of common protein targets (d) Cellular Component of common protein targets (e) Results of the enrichment analysis in a bar chart 18-20 .

Fig. 4 .
Fig. 4. (a) Hexose-6-phosphate dehydrogenase (H6PD) secondary structure.The alpha helix and beta strands seen in the protein structure during the 100-ns simulation study are depicted.The plot shows there 24.16% are alpha helix, 17.94% are beta strands, and totally 42.10% are for secondary structures.(b) Protein secondary structure component (SSE) like as beta strands and alpha helices are constantly seen over the 100 ns simulation.Based on residue index, this graph illustrates how SSEs are distributed across structure of the protein.The SSE formation in every path structure during the MD study is summarized and at the bottom tracks the cumulative SSE assignment for every residue.

Fig. 5 .
Fig. 5. (a) Protein-ligand contact plots (b) The timeline graph displayed the interactions between the Hexose-6-phosphate dehydrogenase with Piperine in all directions in the trajectory frames (a deeper orange color indicates more than one explicit interaction with the ligand).

Table 1 .
Docking Grid Box parameters used in AutoDock Software.

Table 3 .
Minimum binding energy scores of Piperine.