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

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 S128,29.Table 2 Piperine’s physicochemical characteristics.Network pharmacology protein target analysis
Polycystic ovary syndrome − 988 Genes (https://www.disgenet.org/browser/0/1/0/C0032460/) Hyperandrogenism − 108 Genes (https://www.disgenet.org/browser/0/1/0/C0206081/),Oligomenorrhea − 37 Genes (https://www.disgenet.org/browser/0/1/0/C0028949/). We predicted 5 common genes in “Poly Cystic Ovary Syndrome (PCOS)”, “Hyperandrogenism” and “Oligomenorrhea”: NR3C1, PPARG, FOS, CYP17A1, H6PD. Moreover, through intersection analysis of significant pharmacological targets and genes related to PCOS, A collective of five genes has been recognized as to be viable cross-targets for PCOS treatment (Fig. S1b).Gene ontology’s predictionIn 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.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 chart18,19,20.Protein-protein interaction network analysisWe 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 enrichment27. 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 observed30. 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 analysisSeveral 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 PCOS31. Furthermore, Fig. 1e displays the top pathways with statistically significant differences (P < 0.05) as identified by the KEGG enrichment analysis.Molecular docking analysisFollowing 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 Transcription Factor AP-1, Peroxisome Proliferators’ Activated Receptor Gamma (PPARG), Group C Nuclear Receptor Subfamily 1 Member and Member 1 of Family 17 Subfamily A of Cytochrome P450 (CYP17A1)32. The target and piperine’s molecular docking investigation showed that piperine had the greatest binding affinity. The piperine had − 7.96 Kcal/mol (Ki = 1.43 µM), − 8.34 Kcal/mol (Ki = 771.66 nM), -6.42 Kcal/mol (Ki = 19.77 µM), -6.43 Kcal/mol (Ki = 19.34 µM) and − 8.70 Kcal/mol (Ki = 420.17 nM) docking scores respectively show in Table 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.Table 3 Minimum binding energy scores of Piperine.Fig. 2(a) Group C Nuclear Receptor Subfamily 1 Member in Complex with Piperine. (b) Piperine complex with Peroxisome Proliferator Activated Receptor Gamma (c) Transcription Factor AP-1 Subunit C-Fos in complex with Piperine (d) Cytochrome P450 Family 17 Subfamily A Member in complex with Piperine (e) Combined with piperine, hexose-6-phosphate dehydrogenase / glucose 1-dehydrogenase.Fig. 32D Interaction Plot of Docked Complexes (a) Group C Nuclear Receptor Subfamily 1 Member (b) Gamma-Activated Peroxisome Proliferator Receptor (c) Transcription Factor AP-1 Subunit C-Fos (d) Member 1 of Family 17 Subfamily A of Cytochrome P450 (e) Glucose 1-dehydrogenase/hexose-6-phosphate dehydrogenase (H6PD).Molecular dynamics (MD) simulationsThe 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 structure33. 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 molecules34. 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 terms35.The overall potential energy of the system is determined by accounting for all potential interactions between atoms within the protein system.$$\begin{aligned} U\left( {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{R} } \right) = & \sum\limits_{{{\text{bonds}}}} {{\text{K}}_{{\text{b}}} \left( {{\text{b}} – {\text{b}}_{0} } \right)^{2} } + \sum\limits_{{{\text{angles}}}} {{\text{K}}_{{{\theta }}} \left( {{{\theta }} – {{\theta }}_{0} } \right)^{2} } + \sum\limits_{{{\text{dihedrals}}}} {{\text{K}}_{\chi } \left( {\chi – \chi _{0} } \right)} ^{2} \\ & + \sum\limits_{{{\text{nonbond}}}} {\varepsilon _{{{\text{ij}}}} \left[ {\left( {\frac{{\text{R}_{{{\text{ij}}}}^{{\min }} }}{{r_{{{\text{ij}}}} }}} \right)^{{12}} – \left( {\frac{{\text{R}_{{{\text{ij}}}}^{{\min }} }}{{r_{{{\text{ij}}}} }}} \right)^{6} } \right]} + \sum\limits_{{{\text{nonbond}}}} {\frac{{{\text{q}}_{{\text{i}}} {\text{q}}_{{\text{j}}} }}{{{\text{Dr}}_{{{\text{ij}}}} }}} \\ \end{aligned}$$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 dynamics35.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 simulations36. 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 trajectories37. 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 interactions38. 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 metrics39.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 methodologies38.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 resources40. 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 geometries41. 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 simulations41. An RMSD value ≤ 2 Å is considered optimal and recommended42. 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 RMSDThe 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 RMSDThe 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 ligand43.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 absence44. This analysis helps characterize how ligand binding influences the structural dynamics of the protein at the atomic level.Protein RMSF calculationThe 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.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.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:$${{RMSF}}_{{\text{i}}} {\text{ = }}\sqrt {\frac{1}{T}\sum\limits_{{t = 1}}^{T} {\left( {r_{i}^{{\prime}} (t) – r_{i} \left( {t_{{ref}} } \right)} \right)^{2} } }$$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 calculationThe 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.$${{RMSF}}_{{\text{i}}} {\text{ = }}\sqrt {\frac{1}{T}\sum\limits_{{t = 1}}^{T} {\left( {r_{i}^{{\prime}} (t) – r_{i} \left( {t_{{ref}} } \right)} \right)^{2} } }$$where T is the trajectory time over which the RMSF is calculated, tref is the reference time, ri 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 residuesProtein-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.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).Hydrogen bond contactsH-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, side-chain 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.5 Å 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 acceptor45. 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 determined46. 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 hydration47. 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 contactsHydrophobic 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 interactionIonic 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 bridgesWater Bridges are hydrogen-bonded protein-ligand interactions mediated by a water molecule. The hydrogen-bond 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 timelineA 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).