The leaves of Punica granatum L. plant were collected from Gandhi Krishi Vignan Kendra (GKVK), Bengaluru, Karnataka, India. After authentication by the plant taxonomist, a voucher specimen (BOT/mLAC/BMPG001) has been submitted to Botany Department, Maharani Lakshmi Ammanni College for Women, Bengaluru, Karnataka, India. The leaves were shade dried, grounded, and the powder was subjected to soxhlation using 70% methanol as a solvent system. The extract was lyophilized using freeze-drying technology under –55 ${}^{\circ }$C and at constant pressure. The extract thus obtained was used further for analysis.

The yield of hydromethanolic leaf extract was calculated using the following equation: Yield (g/100 g of dry plant material) = (W1 $\times$ 100) / W2, where W1-weight of the extract after solvent evaporation and W2-weight of dry plant material (leaves).

Qualitative phytochemical analysis of HMPGL, for the detection of alkaloids, carbohydrates, flavonoids, glycosides, phenols, saponins, steroids, tannins [19], fatty acids, coumarins and resins [20] was carried out as per the standard procedures.

DPPH radical scavenging assay was used to determine the antioxidant activity of HMPGL. Sample aliquots of different concentrations were incubated with 1.8 mL of DPPH for 30 min at room temperature and the absorbance was measured at 540 nm [19]. The ability of the HMPGL to scavenge DPPH radical was calculated using the standard formula: DPPH free radical scavenging activity (%) = (Absorbance of control (DPPH)–Absorbance of DPPH radical + HMPGL)/(Absorbance of control (DPPH)) $\times$ 100.

The total phenolic content of HMPGL was estimated using the Folin-Ciocalteu reagent as previously reported by Goyal et al. 2010 [19].

The gas chromatogram-mass spectrometer (Agilent Technologies 7890B GC and Triple Quadrupole mass spectrometer 7000D series) was used for the analysis, having a blend silica column, bundled with 5% biphenyl 95% dimethylpolysiloxane (Elite-5MS), merged with a capillary column (30 m $\times$ 0.25 mm $\times$ 1 $\mu$m). Helium quench gas and nitrogen collision gas were deployed as a carrier gas at a regular flow rate of 2.25 mL/min and 1.5 mL/min, respectively, and a pressure of 8.745 psi to separate different components found in the test sample. Each chromatographic run was adjusted at 260 ${}^{\circ }$C injector temperature. 2 $\mu$L volume of the 100$\times$ diluted sample was injected in the instrument (a split ratio of 10 : 1). The oven temperature was in the range of 20 ${}^{\circ }$C to 260 ${}^{\circ }$C. The mass detector conditions were as follows: transfer line temperature 300 ${}^{\circ }$C; ion source temperature 260 ${}^{\circ }$C; ionization mode electron impact at 70 eV, a scan time 0.2 sec and scan interval of 0.1 sec. The sample was analyzed within a mass range of 50 m/z to 1000 m/z. The solvent delay was 0 to 2 min, and the total GC-MS/MS running time was 44 min. The mass spectrum of each peak in the total ion chromatogram was compared with the databases of The National Institute of Standards and Technology (NIST) 14 MS Library containing ‘276248’ number of compound’s spectral databases. The IUPAC names, molecular formulae, molecular weights, and structure of the metabolites of HMPGL were thereby ascertained.

The GC-MS/MS results were obtained to find the most probable lead molecules present in HMPGL. The metabolites thus obtained, if were available in PubChem, then their structures were downloaded from the PubChem compound database (https://www.ncbi.nlm.nih.gov/pccompound) [21] in 3D conformation and .sdf format. The structures of few metabolites, which were not available in PubChem, were obtained from NIST Library and sketched using ChemSketch. The standard protocols of BIOVIA Discovery Studio v3.5 (DS) software (Accelrys Software Inc., USA, 2012) was used for ligand preparation. This standard program was widely employed to prepare all the lead-like compounds to fix and resolve all different chemical properties such as different protonation states, ionization states, isomers, tautomers adding hydrogen, removing duplicates, and fixing bad valencies. This step is crucial because the receptor-ligand interactions have different protonation states; isomers and tautomers typically have different 3D geometries and binding characteristics.

This study is mainly centered on drug target proteins related to breast cancer; therefore, a list of principal biomarkers and prognostic biomarkers that were immunohistochemically tested in breast cancer tissue, are considered for the study [22]. A network of these proteins was constructed using a string database version 11.0 (https://string-db.org/) [23] and analyzed using Cytoscape, and the protein with a medium number of edges was used as a protein target. Basic sequence and structure analysis of the ER was conducted. The protein structures of the different isoforms of ER structure such as ER$\alpha$ LBD (PDB code 3ERD), ER$\beta$ LBD (PDB code 3OLS) and ERR$\gamma$ LBD (PDB code 2GPU) were retrieved in .pdb formats from structural protein database (https://www.rcsb.org/) [24]. Initially, hetatoms and unbound water molecules were removed, following which, a CHARMM force field was applied. It removes alternative conformers and balances valencies of the amino acids. Subsequently, the protein predation protocol builds a loop and refines the side chain of proteins.

The 145 compounds identified through GC-MS/MS analysis were then subjected to a dual-step virtual screening process, in which, first, the compounds were screened for their pharmacokinetic activities, such as ADMET, which were predicted using ADMET Descriptors protocol in DS v3.5.

ADME studies are widely employed in drug discovery to optimize the properties to convert leads into drug molecules or medication. The lead compounds are mostly identified through virtual high-throughput screening approaches or by virtual screening. Drug discovery statistics reveal that around 50% of the drugs fail in the course of clinical trials due to nonstandard effectiveness (efficacy), which has low bioavailability due to poor intestinal absorption and undesirable metabolic stability [25]. The failure of drug-like candidates in the later stages of drug development proves very expensive. Therefore, to scale back the value and clinical failures of recent drug-like molecules, the lead compounds were screened within the initial stages for ADMET. Secondly, the ADMET screened compounds were checked for Lipinski’s rule of five (RO5) violations to ensure that the compounds have drug-likeness properties [5]. Finally, the remaining compounds were taken for further process.

The docking study was performed by identifying the binding site of drug-target protein using the LibDock algorithm available in Discovery Studio 3.5. It uses protein binding site features to direct docking. It finds polar and apolar probes by placing a grid around the ligand 20 Å by 20 Å by 20 Å and extra space of 5 Å in each direction [26]. The site volume and binding site sphere will vary for all the three drug targets and were as follows; for ER$\alpha$ LBD 5.524X-0.284Y-5.750Z, ER$\beta$ LBD 25.17X-27.96Y-11.25Z and ERR$\gamma$ LBD 62.44X-47.17Y-25.727Z. Bad hotspots were removed manually. Finally, pose optimization was done using Broyden-Fletcher-Goldfarb-Shanno (BFGS) and top-scoring ligand poses were ranked and retained.

The molecular dynamic simulation (MDS) was ideally carried out to inspect the stability and rationality of the binding patterns between the probable ligands and the specific target protein. The finest interaction hits obtained from receptor-ligand docking were used for MDS. In this experiment, top pose and elite compound interaction with individual drug targets were simulated using the GROMACS version 2016 package [3, 7]. 3ERD-2-propenoic acid, 3-(4-hydroxyphenyl)-, methyl ester, 3OLS-2-propenoic acid, 3-(4-hydroxyphenyl)-, methyl ester and 2GPU-2-propenoic acid, 3-(4-hydroxyphenyl)-, methyl ester complex topologies were prepared separately for the protein and the ligand. For the ligand molecular topology, the coordinate files were generated by the PRODRG 2.5 server, followed by solvating the receptor-ligand complex in the dodecahedral box with a minimal distance of 1.0 nm. The whole protein-ligand complex and aqueous system are maintained at neutral conditions using Na${}^{+}$ and Cl${}^{-}$ions. Subsequently, the steepest descent algorithm was employed to implement a minimum of 1000 to a maximum of 50,000 steps of minimization to retain the proper distance between the atoms accompanied by no structural changes in the compound. The temperature of 300 K and pressure of 1 bar were attained by making use of canonical ensemble (NVT) and isothermal-isobaric ensemble (NPT) equilibration simulations of 500 ps. The time step of all stages was set to 2 fs. In the end, the protein-ligand complexes were subjected to 30 ns MDS. For further investigation of the stability and flexibility of the complexes, root mean square deviation (RMSD), root mean square fluctuation (RMSF), the radius of gyration (Rg), and hydrogen bonds were analyzed through graphical representation.

The % yield of the HMPGL was found to be 1.44%. Qualitative phytochemical analysis revealed the absence of resins and the presence of alkaloids, carbohydrates, coumarins, flavonoids, glycosides, saponins, steroids, tannins, fatty acids, phenols and terpenoids. The IC${}_{50}$ value for the antioxidant activity of the extract based on the DPPH free radical scavenging assay was identified to be 3.12 mg/L and was comparable with the standard quercetin (3.91 mg/L). Phenol content of the HMPGL found to be 86.17 $\pm$ 3.67 mg gallic acid equivalent 100 g${}^{-1}$ dry weight.

Fig. 1 indicates the GC-MS/MS chromatogram of the HMPGL. After comparing mass spectra of the components with the NIST Library, 145 phytocompounds were identified and listed in Table 1. Out of the 145 compounds screened, 102 and 42 compounds (Table 2, Ref. [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67]) were identified to have more than 80% match with NIST library, respectively. Among the screened compounds, 19 were identified as novel metabolites. The most prevailing compounds among the 145 metabolites are 9,12,15-octadecatrienoic acid, (Z,Z,Z)- (linolenic acid), 1,2,3-benzenetriol(pyrogallol), n-hexadecanoic acid (palmitic acid), maltol, 5-hydroxymethylfurfural, 4H-Pyran-4-one, 2,3-dihydro-3,5-dihydroxy-6-methyl-, undecanol-5 and 3-furaldehyde. These phytocompounds are different fatty acids, terpenes, heterocyclic compounds, flavonoids, pyrrolidines, sesquiterpenoids and phenols. The mass spectrum corresponding to the structure of 2-propenoic acid, 3-(4 hydroxyphenyl)-, methyl ester (also known as 4-coumaric acid methyl ester), which exhibits the highest binding affinity with the breast cancer receptors considered in the current investigation, was represented in Fig. 2.

Fig. 1.

GC-MS/MS chromatogram of the hydromethanolic extract of Punica granatum L. leaves.

Fig. 2.

Mass spectrum of the component (2-propenoic acid, 3-(4-hydroxyphenyl)-, methyl ester) at retention time (RT) 15.2297 (represented in the insert) with a y-axis scale from 0 to 0.9 $\mathrm{\times }$ 10${}^{\mathrm{2}}$ and its structure identified in hydromethanolic extract of Punica granatum L. leaves by GC-MS/MS. Extracted ion chromatogram (EIC) peaks is the same component spectrum at a scale from 0.6 to 2.0 $\times$10${}^6$ (insert). Molecular formula: C${}_{10}$H${}_{10}$O${}_3$, Molecular weight: 178.18.

The protein network of text mined and experimental interactions; and interaction score with 0.9 confidence, consisted of 81 number of nodes and 234 edges with a P-value of 1.0${}^{e-16}$, 0.072, network density and 0.919 network heterogenicity. The network revealed Tumor protein 53 (TP53), RAC-alpha serine/threonine-protein kinase (AKT1), cyclin D1(CCND1), epidermal growth factor receptor (EGFR), phosphatase and tensin homolog deleted on chromosome 10 (PTEN), human epidermal growth factor receptor 2 (ErbB2/HER2), DNA mismatch repair protein (MSH2) and MYC as the hub proteins and ER as a suitable drug target with a medium number of edges, a middle degree node, 0.533 closeness centrality and 0.13 clustering coefficient (Fig. 3).

With the advent of new technology and computational tools, there is a massive increase in the deposition of protein structures in PDB. This consequently creates an increased level of difficulty in selecting an optimal PDB entry for docking. X-ray crystallographic structures of ER$\alpha$ LBD (PDB ID:3ERD), ER$\beta$ LBD (PDB ID:3OLS) and ERR$\gamma$ LBD (PDB ID: 2GPU) were selected after a thorough screening of 329 entries of ER’ structures with a minimum resolution, no missing atoms and completeness of the binding domain. The ER$\alpha$ LBD, ER$\beta$ LBD and ERR$\gamma$ LBD have orthogonal bundle architecture. They mainly consist of $\alpha$-helix, belong to the family of nuclear receptor and exhibited retinoid-X-receptor topology. Largely, there are numerous conserved regions observed in-between the three sequences. However, ER$\alpha$ LBD and ER$\beta$ LBD showed 83.5% sequence and structure similarities, whereas ERR$\gamma$-ER$\alpha$ and ERR$\gamma$-ER$\beta$ showed a similarity of 62.2% and 61.5%, respectively. Overall, multiple sequence alignment showed 55.7% similarity between the three LBDs of ER$\alpha$, ER$\beta$ and ERR$\gamma$. The raw protein structures were refined using a protein preparation tool prior to molecular docking. A loop of ER$\beta$ and ER$\alpha$ was built based on SEQRES data, the atoms were arranged, and its side chains were refined using PDB data. Notably, the structural analysis also revealed that the binding pocket of ER$\alpha$ is larger and broader than that of ER$\beta$ (Fig. 4).

Fig. 3.

Protein interaction network of principal and prognostic biomarkers of breast cancer.

Fig. 4.

Active site representation of isomers of human estrogen receptor.

Efficacy and toxicity are the pivotal determinants of successful drug development. Therefore, all the GC-MS/MS characterized metabolites were subjected to the virtual screening process thoroughly. Initial screening with Absorption, Distribution, Metabolism, Excretion and Toxicity (ADMET) drastically narrowed down the count of molecules (Figs. 5,6) followed by RO5 violation AlogP $\le$5, with molecular weight $\le$500, hydrogen bond donor $\le$5, and hydrogen bond acceptor $\le$10 were evaluated to justify their drug-likeness behaviour. Fig. 5 represents the pharmacokinetic profile of all the metabolites identified in HMPGL. Out of 145 compounds analyzed, 77 revealed good and optimal solubility with a score of 3 and 4. Total, 103 molecules demonstrated high, medium and low blood-brain barrier (BBB) penetration. The human intestinal absorption of 109 compounds was 0 and 1, indicating good and moderate absorption of these hits (Fig. 5). 125 compounds were non-inhibitors of CYPD26, 96 were nontoxic to liver cells and 50 exhibited low binding affinity to plasma proteins. The 16 and 47 compounds showed no carcinogenicity in female and male rats, respectively. 122 compounds were non-mutagens in TOPKAT screening (Fig. 6). Therefore, systematic analysis of all these pharmacokinetic properties revealed only nine compounds that are nontoxic and safe. Among them, 9,12,15-octadecatrienoic acid, 2,3-dihydroxypropyl, ester (Z,Z,Z)- and hexadecanoic acid, 2-hydroxy-1-(hydroxymethyl)ethyl ester were not drug-like molecules as their AlogP were ‘greater than five’ (violation of Lipinski’s RO5). However, 19 novel metabolites were not docked against the ER due to their limitation of exhibited carcinogenicity based on in-silico prediction.

Fig. 5.

The graph of ADMET 2D polar surface area (PSA_2D) vs. ALogP of 145 phytochemicals of Punica granatum L. representing the confidence limit ellipses of 95% and 99% corresponding to the intestinal absorption models and blood-brain barrier (BBB). (a,b) represent molecules before and after ADMET screening. (c) The ellipses define regions where well-absorbed Human Intestinal Absorption (HIA) compounds are expected to be found after oral administration. (d) The graph represents compounds with low and medium blood-brain barrier (BBB) penetration.

Fig. 6.

Number of phytocompounds of hydromethanolic extract of Punica granatum L. leaves predicted to be CYPD26 inhibitors, non hepatotoxic, no plasma protein binding (PPB) affinity, non-carcinogens in Female (F) and Male (M) non-mutagens and obey RO5.

Detailed tabulations of the screened compounds, which completely satisfied the ADMET properties and drug-likeness, are given in Table 3. It is evident from the table that the standard drug, tamoxifen, an anti-estrogen that is widely used in the clinic to treat ER-positive breast tumors) violated both RO5 and exhibits hepatotoxicity and high affinity to plasma binding proteins. 2-Propenoic acid, 3-(4-hydroxyphenyl)-, methyl ester, a natural compound, followed RO5 and did not exhibit any toxicity in ADMET prediction suggesting its efficacy and nontoxic nature. Some of the natural compounds of HMPGL also showed hepatotoxicity and affinity to plasma proteins.

The screened compounds and the standard drug tamoxifen were docked with three different ER structures. The docked poses were examined based on the LibDock score (Kcal/mol) and various types of interactions in hydrogen/hydrophobic interaction analyses. The structural conformations of 4-coumaric acid methyl ester (92203) were the most favorable for the binding cavity of all the three receptors, especially hydroxyl group (-OH) exhibits a major structure-activity relationship, where the removal of (-OH) in the para position of the structure is directly proportional to dock score values. Thus, this position is considered crucial for the binding of this phytochemical to these receptors. Other compounds like benzoic acid, 3,4,5-trihydroxy-, methyl ester (7428), 2-propanone,1-hydroxy-3-(4-hydroxy-3-methoxyphenyl) (586459), (4H)4a,5,6,7,8,8a-hexahydrobenzopyran-5-one-3-carboxamide,2-(2-hydroxypentyl)-8a-methoxy-4a-methyl (587806) and tamoxifen (2733526) also showed favorable binding, but the standard drug tamoxifen violated the RO5 and demonstrated toxicity. The structure of 4-coumaric acid methyl ester showed high complementarity to $\alpha$LBD with a docking score of 75.16 kcal/mol, whereas $\gamma$LBD is favorable for tamoxifen. Nevertheless, $\beta$LBD was structurally compatible with both 4-coumaric acid methyl ester and tamoxifen, as represented in Fig. 7.

Fig. 7.

The binding energy of 4-coumaric acid methyl ester present in hydromethanolic extract of Punica granatum L. leaves and tamoxifen with three isoforms of estrogen receptor (3ERD, 3OLS and 2GPU corresponds to the PDB ID of the LBDs of ER$\alpha$, ER$\beta$ and ERR$\gamma$, respectively).

The binding pattern and chemical interactions of nontoxic natural 4-coumaric acid methyl ester present in HMPGL with good docking score, with the three isoforms of ER, were demonstrated in Fig. 8. The $\alpha$ (W393, E353, R394, L449 and E323) and $\beta$ (F356, L343, L339 and A302) LBD residues are considered as crucial active site amino acids for ligand binding and responsible for antagonist activity. R394 of $\alpha$ binding domain makes crucial Pi interactions with the ligand. In contrast, R316 and L309 amino acids play a vital role in ERR$\gamma$ LBD. In $\alpha$LBD, 4-coumaric acid methyl ester interacts with W393, R353, R394, L449 and E323 to form strong hydrogen interaction and Pi-cation, respectively, whereas, with $\beta$ and related $\gamma$, it forms hydrogen bonding with F356 and R316, and L309, respectively. In all these interactions, -OH forms hydrogen bond interaction; so, it can be considered as an active pharmacophore for 4-coumaric acid methyl ester.

A time-dependent MD simulation at 50 ns was conducted using GROMACS 2016 to investigate the flexibility and overall stability of docked complexes. MD simulation were carried for the best candidate molecule, 4-coumaric acid methyl ester, with three drug target receptors and time-dependent parameters, were analyzed as described by Dhivya et al., 2018 [3]. The root mean square deviation (RMSD), the radius of gyration (Rg) and hydrogen bond (H-bond) interaction graphs were generated by using Xmgrace software in a Linux environment (Fig. 9).

Fig. 8.

Both 3D and 2D docked representation 4-Coumaric acid methyl ester with three different estrogen receptors’ (drug targets) active sites. 4-Coumaric acid methyl ester docked with 3ERD (the ligand-binding domain of ER$\alpha$; A&B), 3OLS (the ligand-binding domain of ER$\beta$; C&D), and 2GPU (the ligand-binding domain of ERR$\gamma$; E&F).

Fig. 9.

Time-dependent MD and simulation of 4-coumaric acid methyl ester and tamoxifen with three different ERs’ (drug targets) active sites (LBD of ER$\alpha$, ER$\beta$ and ERR$\gamma$). (a) Root mean square deviation (RMSD) plot. (b) Radius of gyration (Rg) plot. (c) Hydrogen bond (H-bond) plot.

As is evident in Fig. 9a.1, the increasing trend of RMSD was observed for all the six complexes, 3ERD/ER$\alpha$ LBD–4-coumaricacid methyl ester, 3OLS/ER$\beta$ LBD–4-coumaricacid methyl ester, and 2GPU/ERR$\gamma$ LBD–4-coumaricacid methyl ester, having diverse RMSD values 0–0.30. Initially, the ER$\alpha$ LBD complex showed little jiggling, whereas ER$\beta$ LBD showed high fluctuations, but ERR$\gamma$ LBD exhibited a gradual increase in RMSD until 10 ns. The fluctuations of the ER$\alpha$ LBD complex were between 0.10 to 0.30 nm, but after 30 ns, this complex attained stability with slight deviations in RMSD value. However, the fluctuations of the ER$\beta$ LBD complex elevated up to 0.25 nm at 5 ns, dropped to 0.2 nm at 12.5 and 20 ns, with a sudden drop of RMSD to 0.2 nm at around 12.5 and 20 ns, fluctuated till 35 ns and steadily increased at the end of dynamics. The RMSD of the ERR$\gamma$ LBD complex gradually increased until 35 ns and attained plateau, which was in the range of 0–0.2. The Rg analysis usually measures the compactness of a protein or system. The predicted results in Fig. 9b.1 demonstrated that all the three complexes ER$\alpha$ LBD, ER$\beta$ LBD and ERR$\gamma$ LBD- were static and constant with Rg values 1.8 to 1.85 nm, 1.775 to 1.825 and 1.750 to 1.850 nm, respectively, throughout the simulation time 0–50 ns, the complexes slightly fluctuated to attain its stability. Comparative MD simulations revealed that the residual backbone and appropriate conformation of the ER$\alpha$ LBD complex was stable compared to other complexes. The H-bond plot of the simulation was presented in Fig. 9c.1. It is practically impossible to inspect H-bond stability by examining the crystal, nuclear magnetic resonance (NMR), and electron microscopy-solved structures. Consequently, it is significant to probe the stability of the H-bonds utilizing MD simulations. The presence of H-bond interactions in the docked complexes was identified by the gmx hbond tool in the accepted geometry (distance 3.5 Å and 180${}^{\circ }$$\pm$ 30${}^{\circ }$, between the donor and acceptor). H-bonds between the ligand and the receptor jiggled all through the simulation. The maximum number of H-bonds bound confirmation for ER$\alpha$ LBD, ER$\beta$ LBD and ERR$\gamma$ LBD complexes are 3, 3 and 4, respectively. It is evident from the 3ERD/ER$\alpha$ LBD–Tamoxifen, 3OLS/ER$\beta$ LBD–Tamoxifen, and 2GPU/ERR$\gamma$ LBD–Tamoxifen, having diverse RMSD values 0–0.30. (Fig. 9a.2). All the three complexes exhibited higher stability compared to complexes represented in Fig. 9a.1. From Fig. 9b.2 it can be noted that 3OLS/ER$\beta$ LBD–Tamoxifen was more stable all through the run. Only 3OLS/ER$\beta$ LBD–Tamoxifen showed H interactions.