Gas chromatography-mass spectrometry metabolic profiling, molecular simulation and dynamics of diverse phytochemicals of Punica granatum L. leaves against estrogen receptor

Introduction: Breast cancer is the most common type of cancer globally and its treatment with many FDA-approved synthetic drugs manifests various side effects. Alternatively, phytochemicals are natural reserves of novel drugs for cancer therapy. Punica granatum commonly


Abstract
Introduction: Breast cancer is the most common type of cancer globally and its treatment with many FDA-approved synthetic drugs manifests various side effects. Alternatively, phytochemicals are natural reserves of novel drugs for cancer therapy. Punica granatum commonly known as pomegranate is a rich source of phytopharmaceuticals. Methods: The phytoconstituents of Punica granatum leaves were profiled using GC-MS/MS in the present work. Cytoscape-assisted network pharmacology of principal and prognostic biomarkers, which are immunohistochemically tested in breast cancer tissue, was carried out for the identification of protein target. Followed by, rigorous virtual screening of 145 phytoconstituents against the three ER isoforms (α, β and γ) was performed using Discovery Studio. The docked complexes were further evaluated for their flexibility and stability using GRO-MACS2016 through 50 ns long molecular dynamic simulations. Results: In the current study, we report the precise and systematic GC-MS/MS profiling of phytoconstituents (19 novel metabolites out of 145) of hydromethanolic extract of Punica granatum L. (pomegranate) leaves. These phytocompounds are various types of fatty acids, terpenes, heterocyclic compounds and flavonoids. 4-coumaric acid methyl ester was identified as the best inhibitor of ER isoforms with drug-likeness and no toxicity from ADMET screening. γ-ligand binding domain complex showed the best interactions with minimum RMSD, constant Rg, and the maximum number of hydrogen bonds. Conclusion: We conclude that 4-coumaric acid methyl ester exhibits favourable drug-like properties comparable to tamoxifen, an FDA-approved breast cancer drug and can be tested further in preclinical studies.

Introduction
Drug discovery is a multifaceted and interdisciplinary endeavor that pursues a sequential process, begins with the discovery of a target and lead, followed by lead optimization and preclinical studies to define and verify the suitability of such lead agents through a number of predetermined guidelines for kick-starting clinical development [1,2]. The high cost and time-consuming nature of these processes against an ailment demand a novel kind of costeffective and less time-consuming, intensive in-silico approach [2,3]. In-silico techniques employ bioinformatics tools to identify drug targets, to explore target protein structure for prediction of possible active sites to dock molecules with targets, rank ligands based on their binding affinities (i.e., energies), and optimize lead molecules [4]. Proficiency in calculating precise, significant, and diverse conformations of ligands within the target site make the molecular docking the most desired drug designing tool [5]. Hence, this tool can be used to understand the various interactions between phytochemicals and breast cancer targets or biomarkers. Despite the availability of modern tools and increasing expertise to identify a new anticancer molecule or lead, the potential of natural products (medicinal plants) and their compounds remains immense and undiscovered.
Many FDA-approved drugs such as paclitaxel and vinblastine (anticancer drugs), and quinine and artemisinin (antimalarial drugs) are plant-based. Therefore, combining the potential of phytomedicine with modern technology can help to develop more effective and economical drugs [6]. Phytochemicals from medicinal plants not only serve as drugs but also provide a lead for the development of new drugs. Punica granatum L., commonly known as pomegranate, is one such medicinal plant with a known reservoir of active biomolecules such as a series of compounds known as ellagitannins, and specifically punicalagin, the unique largest molecular weight polyphenol known to date [7]. A plethora of research reports describe the vasculo protective [8], antidiabetic, antioxidant [7], antiproliferative [9], antitumor [9,10], anti-inflammatory [11], neuroprotective [12], and antimicrobial [13] properties of pomegranate along with its potential for treating asthma [14]. Previously, Usha et al. (2014) [15] have identified anticancer targets from pomegranate using an in-silico approach.
Breast cancer is the most common type of cancer worldwide, which affects 2.1 million women annually.
In the year 2018, the World Health Or-ganization estimated that 627,000 deaths in women were caused due to breast cancer accounting for nearly 15% of all cancer-related mortality among women (https://www.who.int/cancer/detection/breastcancer/en/) [16]. The current study aims at identifying a lead molecule for breast cancer from a natural source like pomegranate.
Aromatase is an enzyme that converts C19 androgen to C18 estrogen that plays a vital role in breast cancer, which was reported to be present in high concentrations in breast tumors [17]. Several ellagitannins-derived compounds isolated from pomegranate were previously reported to have potential anti-aromatase activity and inhibit testosterone-induced breast cancer cell proliferation [18]. However, there is a dearth of knowledge about the molecular basis and mechanism(s) of these compounds against breast cancer. The aim of the present study is to design a drug-like entity for breast cancer from Punica granatum L. phytochemicals. To be specific, our study involves a gas chromatography-mass spectrometry (GC-MS)/MS-based metabolite profiling of hydromethanolic extract of Punica granatum L. leaves (HMPGL), to identify a lead candidate, which can act as an antagonist and further verification against the ligand-binding domains (LBDs) of the three different isoforms of human estrogen receptor (ER), namely, ER alpha (ERα), ER beta (ERβ) and estrogen-related receptor (ERR) gamma (ERRγ) using molecular dynamics approach. Our finding reveals that 2-propenoic acid, 3-(4-hydroxyphenyl)-, methyl ester alias 4-coumaric acid methyl ester could act as an antagonist for LBDs of ERα, ERβ and ERRγ isoforms of ER, suggesting its future potential for the development of breast cancer therapeutics.

Collection of plants and preparation of crude extract
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 freezedrying technology under -55 • C and at constant pressure. The extract thus obtained was used further for analysis.

Determination of hydromethanolic leaf extract yield
The yield of hydromethanolic leaf extract was calculated using the following equation: Yield (g/100 g of dry plant material) = (W1 × 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.

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

Phytoconstituent profiling of hydromethonolic extract by GC-MS/MS analysis
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 × 0.25 mm × 1 µ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 • C injector temperature. 2 µL volume of the 100× diluted sample was injected in the instrument (a split ratio of 10 : 1). The oven temperature was in the range of 20 • C to 260 • C. The mass detector conditions were as follows: transfer line temperature 300 • C; ion source temperature 260 • 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.

Preparation of ligands and proteins
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α LBD (PDB code 3ERD), ERβ LBD (PDB code 3OLS) and ERRγ 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.

Absorption, Distribution, Metabolism, Excretion and Toxicity (ADMET) and drug-likeness screening
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 AD-MET, 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 AD-MET. 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.

Protein-ligand docking
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α LBD 5.524X-0.284Y-5.750Z, ERβ LBD 25.17X-27.96Y-11.25Z and ERRγ 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.

Molecular dynamics and simulation (MDS)
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-2propenoic 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.

Sequence and structure assessment of estrogen receptors and their related ligand-binding domains
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 Pvalue 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α LBD (PDB ID:3ERD), ERβ LBD (PDB ID:3OLS) and ERRγ 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α LBD, ERβ LBD and ERRγ LBD have orthogonal bundle architecture. They mainly consist of α-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α LBD and ERβ LBD showed 83.5% sequence and structure similarities, whereas ERRγ-ERα and ERRγ-ERβ showed a similarity of 62.2% and 61.5%, respectively. Overall, multiple sequence alignment showed 55.7% similarity between the three LBDs of ERα, ERβ and ERRγ. The raw protein structures were refined using a protein preparation tool prior to molecular docking. A loop of ERβ and ERα was built based on SE-QRES 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α is larger and broader than that of ERβ (Fig. 4).

In-silico pharmacokinetics properties and Lipinski's rule (RO5) validation
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 ≤5, with molecular weight ≤500, hydrogen bond donor ≤5, and hydrogen bond acceptor ≤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 bloodbrain 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,3dihydroxypropyl, ester (Z,Z,Z)-and hexadecanoic acid, 2hydroxy-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.
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 ta-ble 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 AD-MET prediction suggesting its efficacy and nontoxic nature. Some of the natural compounds of HMPGL also showed hepatotoxicity and affinity to plasma proteins.

Binding and interaction analysis
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-5one-3-carboxamide,2-(2-hydroxypentyl)-8a-methoxy-4amethyl (587806) and tamoxifen (2733526) also showed favorable binding, but the standard drug tamoxifen violated the RO5 and demonstrated toxicity. The structure of 4coumaric acid methyl ester showed high complementarity to αLBD with a docking score of 75.16 kcal/mol, whereas γLBD is favorable for tamoxifen. Nevertheless, βLBD was structurally compatible with both 4-coumaric acid methyl ester and tamoxifen, as represented in Fig. 7. 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 α (W393, E353, R394, L449 and E323) and β (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 α binding domain makes crucial Pi interactions with the ligand. In contrast, R316 and L309 amino acids play a vital role in ERRγ LBD. In α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 β and related γ, it forms hydrogen bonding with F356 and R316, and L309, respectively. In all these interactions, -OH forms hydrogen bond interac-tion; so, it can be considered as an active pharmacophore for 4-coumaric acid methyl ester.

Time-dependent parameter analysis
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, 4coumaric 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).     As is evident in Fig. 9a.1, the increasing trend of RMSD was observed for all the six complexes, 3ERD/ERα LBD-4-coumaricacid methyl ester, 3OLS/ERβ LBD-4coumaricacid methyl ester, and 2GPU/ERRγ LBD-4coumaricacid methyl ester, having diverse RMSD values 0-0.30. Initially, the ERα LBD complex showed little jiggling, whereas ERβ LBD showed high fluctuations, but ERRγ LBD exhibited a gradual increase in RMSD until 10 ns. The fluctuations of the ERα 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β 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γ 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α LBD, ERβ LBD and ERRγ 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α 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 Hbond interactions in the docked complexes was identified by the gmx hbond tool in the accepted geometry (distance <3.5 Å and 180 • ± 30 • , 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α LBD, ERβ LBD and ERRγ LBD complexes are 3, 3 and 4, respectively. It is evident from the 3ERD/ERα LBD-Tamoxifen, 3OLS/ERβ LBD-Tamoxifen, and 2GPU/ERRγ 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β LBD-Tamoxifen was more stable all through the run. Only 3OLS/ERβ LBD-Tamoxifen showed H interactions.

Discussion
Breast cancer is regarded as one of the major burdens observed in women. Consequently, there is an upsurge in efforts for the discovery of new drugs to combat this disease. Plants are widely regarded as a reservoir of various types of bioactive metabolites with various therapeutic and pharmacological potentials [68]. Due to their diverse nature, several plant-based molecules are used as drugs and are in the process of discovery routes [69]. Most of the research studies reported to date involved the characterization of peel, seed and bark by GC-MS and very few studies on leaves of Punica granatum L. using the NMR technique. However, no studies have been reported on leaves by GCMS analysis. Notably, in the present study, we attempted to characterize the possible phytochemicals exclu-sively by GCMS in Punica granatum L. leaves [70,71]. The present study identified 145 phytoconstituents from HMPGL by GC-MS profiling. Similar to our report, 5hydroxymethylfurfural was reported in high concentrations in ethyl acetate extract of the fruit peel of Punica granatum by Barathikannan et al. (2016) [71]. These phytochemicals are known to contribute to the diverse medicinal properties of the plant, as previously reported [69]. A majority of these phytocompounds were reported to have antioxidant, antiinflammatory, anticancer, and antimicro-bial activities, and used as food additives (Table 1). Due to the presence of alkaloids, coumarins, flavonoids, glycosides, phenols, saponins, tannins and terpenoids, HMPGL has better DPPH radical scavenging activity. Bekir et al. (2013) [72] had previously reported the antioxidant potential of Punica granatum L. leaves extracted with different solvents based on their polarity.
Using the in-silico approach, the pharmacokinetics and toxicity studies of compounds are investigated before evaluating their biological activity [73,74]. Poor pharmacokinetic profile and toxicity are the main reasons for last stage failures in drug discovery. Therefore, in the present study, all the identified 145 phytochemicals from the HMPGL were initially screened for ADMET and Lipinski's RO5, and only eight molecules were identified to be drug-like nontoxic molecules. Among the 138 phytoconstituents, a majority of them were identified to be fatty acids and terpenes. Due to the presence of the long hydrophobic -acyl chains and terpenes, the fatty acids exhibited low solubility and high affinity to plasma binding proteins indicating low efficacy. They will be highly toxic as they can penetrate BBB, bind to CYPD26 and exhibit hepatotoxicity.
Docking studies were carried out with 51 metabolites (eight drug-like nontoxic molecules + 42 mentioned in Table 1). Those 42 metabolites were also used for docking studies as most of them exhibited hepatotoxicity, and it is one very common toxic property noted in many FDAapproved drugs like tamoxifen. One of the most important steps in developing a new lead molecule is target selection. Systematic analysis of protein interaction networks (disease networks) and identification of hub protein enhances the understanding of the molecular basis of the disease. They also help in determining the key node as a potential target protein in the drug discovery process. Several hub proteins have been identified to be not suitable as drug targets as their inhibition may affect crucial activities of the cell; However, ER, a middle degree node and with high cluster coefficient, is an optimal drug target for breast cancer treatment [75,76]. Our binding pocket analysis revealed that ERα is larger and broader than ERβ, which has been previously speculated to be due to the possible sequence diversity of ER Paterni et al. (2014) [77]. Targeting more than one receptor helps in an in-depth understanding of the efficacy and toxicity of the drugs, including complex interactions [78]. So, in this study, we exclusively focused on three significant receptors involved in breast cancer mechanisms (ERα, ERβ and ERRγ). In our docking results, there was notable variation seen in the binding energies with the three receptors for all the docked molecules and these differences can be attributed to selective binding of the phytoconstituents to the diverse LBDs of these receptors. In all the interactions, the amino acids Arg, Leu and Glu of the LBD were making crucial interactions with the ligand (4-coumaric acid methyl ester). Interestingly Arg394 alpha binding domain has been reported to make crucial pi interactions with the ligand [79]. The deletion of the -OH of the 4-coumaric acid methyl ester reduces the binding affinity. Hence, it can be considered a pharmacophoric feature due to its hydrogen bonding with the receptors.
MD simulations of the docked complexes from docking studies help refine docking and enhance the accuracy of the binding affinity predictions [3,68]. Post docking MD simulations at 50ns reflected the time-dependent behaviour of the docked complexes. The ligand in the ERRγ LBD complex was having the least RMSD values and attained stability at 30 ns and remained stable till the end compared to ERβ LBD and ERα LBD complexes. The Rg values are almost in the same range (1.75 to 1.85) for all the complexes, whereas the ERRγ LBD complex with minimal vibrations, fluctuations, maximum H-bonds depicts its greater stability compactness and binding affinity with 4coumaric acid methylester. This compound is a natural phenylpropenoic acid compound, belonging to a class of phenolic acids and is a phenolic antioxidant. The antioxidant potential of this compound might contribute to the anticancer effect against breast cancer plausibly due to its high binding affinity to ERα, ERβ and ERRγ LBDs.

Conclusions
To conclude, the current research work was an attempt to draw insight into the structural, functional, and dynamical aspects of phytocompounds of Punica granatum L. In the process, the GC-MS/MS metabolite profiling revealed the presence of 145 compounds. These can be deployed to discover novel drugs against various cancers, as Punica granatum L. is reported to have anticancer properties. Estrogen receptor plays a crucial role in cellular proliferation and differentiation of breast cancer cells, which was identified as a novel target for breast cancer by network pharmacology. 96% of the phytoconstituents exhibited toxicity in the virtual ADMET screen, and 35% did not exhibit drug-likeness. In-silico, molecular docking was performed against three isoforms of ER and compared with the standard drug tamoxifen. 4-Coumaric acid methyl ester, a nontoxic natural, demonstrated the highest affinity with core residues of 3ERD (ERα LBD), with 3OLS (ERβ LBD), with 2GPU (ERRγ LBD). Thus, the docking poses of the identified hit, 4-coumaric acid methyl ester, were further evaluated by dynamic simulation at 50 ns. Among the various ER receptors, the 2GPU complex, which is ERRγ LBD, is "the best" with minimum RMSD, constant Rg and maximum number of H-bonds indicating a stable and compact system. The presence of 4-coumaric acid methyl ester and various other bioactive metabolites reported justifies the use of Punica granatum L. for treating breast cancer by traditional practitioners. Taken together, our results represent a promising starting point for the preclinical evaluation of 4coumaric acid methyl ester as a possible treatment for breast cancer.

Author contributions
SKM and KRS conceived the idea. TU performed in-vitro and computational studies and drafted the manuscript. DS assisted in in-silico experimentations. AKG helped in the experimental procedure and partially drafted the manuscript. HSY helped in the computational facility and wrote a part of the manuscript. DB helped in GC-MS/MS analysis and thoroughly revised the manuscript. SKM arranged the funds and supervised the whole study, edited, and upgraded the final version of the manuscript. All the authors analyzed, discussed the results, and approved the version to be submitted.

Ethics approval and consent to participate
Not applicable.

Acknowledgment
The authors acknowledge Neoscience Labs Private Limited, Chennai, India for allowing to utilize their GC-MS facility and their technical support. DBT-BIF computational facility and BiSEP facility at MLACW was used to carry out the research work. Dhivya also expresses her sincere thanks for the fellowship provided by the DBT-BIF facility at MLACW by Govt of India. The authors are grateful to Dr. V.R. Devraj, Professor, Bangalore Central University and Dr. C.S. Karigar, Professor, Bangalore University, India, for their valuable suggestions. This publication was supported by the Deanship of Scientific Research at Prince Sattam Bin Abdulaziz University, Al-kharj, Saudi Arabia.