Hibiscus anthocyanins can selectively modulate estrogen receptor activity with favorable toxicology: A computational analysis. Yahyea

The estrogen hormone receptor (ER) mediated gene expression mainly regulate the development and physiology of primary and secondary reproductive system alongside bone-forming, metabolism and behaviour. Over-expressed ER is associated with several pathological conditions and play a key role in breast cancer occurrence, progression and metastasis. Hibiscus sabdariffa L. or roselle is a rich source of naturally occurring polyphenolic compounds including anthocyanins and reportedly have strong estrogenic activity. To validate these findings, we have investigated the estrogen receptor binding affinity and safety of some previously recorded polyphenols using a suite of computational methods. Our investigation showed the estrogen-receptor binding potential of Pelargonidin, Delphinidin, Cyanidin, and Hibiscetin are more efficient than popular breast cancer drugs, Tamoxifen and Raloxifene, with favourable toxicological parameters and low half maximal inhibitory concentration. This is the first report to investigate the phytochemical basis of estrogenic activity of Hibiscus sabdariffa L.

Several small-molecule inhibitors or modulators of estrogen receptor known as selective estrogen receptor modulators (SERM) are currently used in the treatment of breast cancer [8] [9]. Popular breast cancer chemopreventive agents like tamoxifen, raloxifene and toremifene are selective estrogen receptor modulator that binds to ER and acts as agonists or antagonists of estrogen to inhibit its proliferative actions on the mammary epithelium [8] [9]. Although this is a favorable strategy to halt cancer progression however, the toxicity associated with the synthetic molecules limits their use in chemotherapy [10] [11]. Thus, plantderived ER agonists or antagonists received considerable attention because of their estrogen-like activity and selective low toxicity [12] [13].
The plant, Hibiscus sabdariffa L. (HS) or roselle is a rich source of naturally-occurring phenolic compounds like alkaloids, flavonoids, anthocyanins, coumarins and phytosterols [14][15] [16] [17]. The alcoholic extracts of the plant showed mild to moderate estrogenic effects on the uteri of immature female rats and exhibited significant effects on the circulating levels Follicle-stimulating hormone (FSH), Prolactin, Estradiol and Testosterone in male Wistar rats [18] [19] [20]. Previous studies reported a strong inhibitory effect of its extracts on the growth and survival of breast cancer cells (MCF-7) and prostate cancer cells (LNCaP) however, the mechanism of inhibition was not well-depicted [21][22]. Although HSextracts reportedly alters the cellular localization of the estrogen receptor, no significant progress has been made to understand the role of its phytochemicals underlying these activities [23]. Previously, we have recorded all the phytochemicals extracted from Hibiscus sabdariffa L. and reviewed their implication in breast cancer chemotherapy [24]. Based on our previous observation, we have selected twelve prominent phenolic compounds from Hibiscus sabdariffa L. (HSP) that could be a potential modulator of estrogen receptor activity. We approached computational techniques as it often used to complement experimental findings that help to prioritize the class of compounds to screen and dramatically reduces the quantity of compounds to screen. The experimental compounds include the five major anthocyanins (Delphinidin, Cyanidin, Malvidin, Peonidin, Pelargonidin), five phytosterols (β-sitosterol, Stigmasterol, α-Spinasterol, Campesterol, Ergosterol) and two other phenolic compounds (Protocatechuic acid and Hibiscetin) which were selected to evaluate their affinity and specificity towards binding and/or modulating estrogen receptor activity using computational docking analysis. Based on the result obtained from the protein-ligand docking studies, toxicology and drug-likeliness of the best-docked compounds were analysed and the 2D quantitative structure-activity relation (QSAR) analysis was performed for these compounds to predict their biological activity.

Preparation of Experimental Receptors (ERα) and Ligands (HSP):
The validated X-ray diffraction structure of human estrogen receptor alpha (ERα) was retrieved from the RCSB Protein Databank (www.rcsb.org) in .pdb format having PDB_ID: 1X7R (2.00 Å resolution). This particular protein structure was selected as it is already co-crystallized with one of the well-known plantbased ER-agonist, Genistein. The 2D/3D structures and the SMILE strings of twelve experimental molecules were retrieved from the NCBI PubChem (www.pubchem.ncbi.nlm.nih.gov) database [25]. The SMILE strings were converted to .MOL files using Open Babel 3.1.1, a freeware used for conversion of chemical formats [26]. Alongside, two already known SERM molecules: tamoxifen and raloxifene, and the biological ligand of ER, 17β-estradiol were considered as the control for this study. The stereochemical clashes of the compounds were removed by energy minimization to get optimized bond distance, bond angles and set dihedrals using Merck molecular force field (MMFF) in VLIfeMDS 4.6 software. The 3D structure of the target protein (ERα) was optimized using YASARA Energy Minimization Server (http://www.yasara.org/minimizationserver.htm) [27].

Molecular Docking Analysis:
The receptor-ligand molecular docking is a dynamic tool for visualizing and understanding the interaction between a receptor protein and its ligand. The molecular docking analysis was performed using the licensed version of the FlexX™ v2.1.6 included in LeadIT™ package software [28]. The active pocket was defined using the binding site of the known cocrystallized reference ligand by including the amino acid residues within a radius of 6.50Å using the protein preparation wizard of the FlexX LeadIT. Further, the reference ligand and crystal water molecules in the binding pocket were removed from the target protein, as evidently they did not play any role in the ligand-binding ability of ERα. The experimental ligands and the control ligands were then docked with the target proteins which is followed by an analysis of hydrogen bonding parameters. FlexX primarily uses incremental construction algorithm that considers the physicochemical features and variable flexibility of compounds while fitting it in the active site of the rigid target protein [29].
To cross-validate the FlexX results, re-docking was performed using a different docking algorithm in the SwissDock web-server (http://www.swissdock.ch/) that uses CHARMM-based software EADock-DSS [30]. The strength of the predicted crystallographic binding orientations is expressed in term of the energy function, where the lowest energy score is considered as the most stable conformation of the proteinligand interaction [29]. For re-docking, energy-minimized ERα protein structure was loaded in the swissdock server and the ligands were selected from the ZINC database [31]. The SwissDock generates almost 500 docking poses for the protein-ligand complex, among which best-docked pose with lowest energy score was selected from the ViewDock plugin of UCSF Chimera [32].

Prediction of ADME/Tox profile & Bioavailability:
The ADME (Adsorption, Distribution, Metabolism and Excretion) and toxicity profiles of the best four hits were performed using PreADMET (https://preadmet.bmdrc.kr), a web server from Yonsei University, Republic of Korea [33]. The PreADMET server uses in vitro results and considers several parameters for evaluating the intestinal absorption of potential drug candidates. Among them, Caco2 (human colon adenocarcinoma cells) model and MDCK (Madin-Darby canine kidney) cell model has been recommended as a reliable in vitro model for the prediction of oral drug absorption [34]. Besides, HIA (human intestinal absorption) data and skin permeability data can be useful in predicting and selecting between oral and/or transdermal drug delivery [35]. Moreover, skin permeability is a risk assessment parameter for all chemicals that may accidentally come into contact with skin [35]. The Plasma Protein Binding (PPB) data and Blood-Brain Barrier Penetration (BBB) data is useful to predict the availability of drug molecule to interact with its pharmacological target after diffusion or transport across cell membranes and the probability of interfering with the central nervous system [36]. PreADMET server uses a genetic functional approximation to select relevant 2D descriptors, followed by Resilient back-propagation (Rprop) neural network analysis to predict ADME data.
Besides several small molecules fails to attain a good drug-like property after confirming a promising ADME/Tox profile because of low bioavailability and poor pharmacokinetics. Thus, pharmacokinetic parameters like gastrointestinal absorption and brain access play a key role in determining bioavailability and brain-toxicity of molecules under investigation. The Brain or IntestinaL EstimateD permeation method (BOILED-Egg) is one such powerful predictive model developed by Daina A et al. [37] in 2016 which counts on the lipophilicity (n-octanol/water partition coefficient or WLOGP) and polarity (Topological Polar Surface Area or TPSA) of small molecules to accurately predict their passive gastrointestinal absorption and brain penetration. The BOILED-Egg predictive analysis was performed using SwissADME server (http://www.swissadme.ch/) [38]of the Swiss Institute of Bioinformatics. The required physicochemical descriptors were generated by providing the SMILE strings of the experimental molecules that translated into a graphical output of the model.

2D Quantitative Structure-Activity Relation (QSAR) studies:
Statistically-based Quantitative Structure-Activity Relation (QSAR) study is a powerful computational tool used in medicinal chemistry that helps to predict the biological activity of chemical compounds based on their structural features. For QSAR studies, 50 already known inhibitors of Er-α subtype and their corresponding inhibitory concentrations (IC50) values (Supplementary, Table 2) were collected from BindingDB database (www.bindingdb.org) [39] managed by the University of California, USA.
Compounds that showed ERα inhibition/modulation in cell-based assays involving ER+ human ductal breast cancer cells (MCF7) with plant product-inspired scaffolds and two to five fused ring skeletons were considered for QSAR model generation and validation. The IC50 values were converted to their respective negative logarithmic values to reduce their numerical difference, which were further used as the dependent variable for QSAR model development. The whole dataset was divided randomly into 'training set' molecules for model generation and 'test set' molecules for model validation (compounds marked in Supplementary Table 2). The random division of training set and test set molecules was based on a trial and error method that generated 50 2D-QSAR models and the process was repeated until a satisfactory r 2 value was attained for statistical significance. The 2D-QSAR analysis was performed using VLifeMDS 4.6 [40] software suite that has multiple variable selection options and regressions methods developed by NovaLead Pharma. The builder module of the VLifeMDS drew the 3D structures of the 50 known molecules which were then energy minimized using the Merck molecular force field (MMFF). To create a 2D-QSAR model, 705 physicochemical and topological descriptors of the molecules were generated. After discarding invariable descriptors (values that are constant for all compounds), 125 physicochemical and 240 topological descriptors were considered. The multiple regression analysis was performed considering Log10(IC50) as the dependable variable and the molecular descriptors as independent variables in a stepwise forward-backward variable selection protocol. Following the model building, the fitness plot for training and test set molecules along with the contribution of the individual descriptor for activity was analysed.
The bioactivity of the experimental molecules was predicted based on the QSAR model using the Generic Prediction module of VLifeMDS.

Data set:
The target protein, estrogen receptor alpha (ERα) is one of the major therapeutic targets for breast cancer. The three-dimensional structure of the protein availed from the protein data bank and was energyminimized using YASARA server. The initial energy of 630912 kJMol -1 of ERα structure was optimized to -151332 kJMol -1 , which enhanced the post-docking results and the reliability of the LeadIT docking protocol. The stereochemical clashes of the twelve experimental ligands ( Table 1) were minimized using the Merck molecular force field (MMFF) to get validated results especially in further studies like molecular docking simulation, toxicology prediction and 2D-QSAR analysis.

Molecular Docking Analysis:
Recently, a few reports[41] [42][43] on the discovery of novel ERα antagonist/agonist using FlexX docking program have been published that reflects the efficiency of FlexX docking algorithm in studying the hormone-receptor interaction. The molecular docking simulation revealed that most of the experimental ligands exhibiting strong affinity towards the core cavity of ERα active site like estradiol with functional residues. In the top four evaluations, three anthocyanins (Pelargonidin, Delphinidin, Cyanidin) and Hibiscetin exhibited the most stable conformation (Figure 2) with ERα with ΔG binding energies -25.9934, -24.2893, -23.5567, -20.7123 kcal.mol −1 respectively. The energy values of these compounds are much lower than the control ligands Tamoxifen and Raloxifene (-16.4166 and -21.0154 kcal/mol −1 respectively), suggesting energy favourable bonding and higher stability. All these compounds framed three or more hydrogen bonds with the core cavity residues of ERα active site alongside the several weak hydrophobic interactions (Van der Waals forces). The total binding affinity score is an energy function which is expressed as a sum of the energy contributed by all matched interacting groups, the lipophilic contact area, ambiguous contact area, clash penalty, and the conformational entropy score of the ligand. The total ERαbinding energy scores of the best-matched poses and the contribution of different energy parameters are given in Table 2. A comparative 2D interaction map of best-matched compounds with the key active site residues of ERα is provided in the Figure 3. Among the anthocyanin ligands, Pelargonidin showed the most stable pose, forming three hydrogen bonds with active site residues Glu353, Leu387 and His524, with bond energies ranging from -4.7 to -3.9 Kcal.mol -1 . Besides, Pelargonidin also exhibited several weak Van der Waal interactions with active site residues. All the anthocyanins exhibited strong binding affinities towards ERα, which was more favourable than the control ligands in term of the energy values. A detailed hydrogen-bonding pattern and weak interactions of other best-matched anthocyanins are provided in Table   3. Among the flavonoid compounds, Hibiscetin expressed the most promising results by forming five Hbonds with Pro324, Glu353, Ile386 and, Trp393 residues of ERα active site, with bond energies ranging from -4.7 to -3.2 Kcal.mol -1 . The total energy score of Hibiscetin was higher than Raloxifen but lower than Tamoxifen, thus Hibiscetin was considered for further analysis. Although ERα is a steroid hormone receptor, surprisingly the five studied phytosterols (Campesterol, Ergosterol, Stigmasterol, α-Spinasterol, β-Sitosterol) did not show favourable interaction with ERα (Supplementary, Table 1). In the validatory re-docking step using SwissDock, similar results (-8.60, -8.58 and -8.8 Kcal.mol -1 respectively for pelargonidin, delphinidin and cyanidin) were observed as the anthocyanins exhibited much lower binding energy values than control ligands (-7.93 Kcal.mol -1 for Tamoxifen and -7.54 Kcal.mol -1 for Raloxifen)( Table 3). The 2D structures of the best matched compounds are shown in Figure 1.

ADME/Tox Profile and Bioavailability:
According to the ADME prediction of the best hits, Pelargonidin and Cyanidin were predicted to be wellabsorbed by Human intestinal cavity with HIA values 83.58% and 72.5% (Table 4). However, Delphinidin and Hibiscetin showed moderate to low intestinal absorbance with HIA ranged between 20-65%. The Blood-Brain Barrier(BBB) Penetration of all four molecules was predicted to be < 1, which signifies CNS (Central nervous system) inactivity. These findings were again validated in BOILED-Egg Model for Gastrointestinal absorption and Brain penetration, which ensured similar results (Figure 4) Table 4). The rate of skin permeability of the molecules expressed in term of Human skin permeability coefficient (LogKp) whose value generally ranges between -3 to +6 for good absorption via the skin, showed low skin permeability of all the studied molecules.

Quantitative Structure-Activity Relation (2D-QSAR) Model & Bioactivity prediction:
In the present study, 50 2D-QSAR models were generated for Log10(IC50) or Bioactivity of 50 already known ERα agonist/antagonist compounds involving total 365 physicochemical and alignment independent topological descriptors by randomization of training and test set data. Among these models, one with the Recently; the toxicology prediction server PreADMET was used in several published reports [46][47] [48] that suggests the reliability of the server. Beside, these in silico toxicology and bioavaibility analysis methods are now widely accepted by pharmaceutical industries to support various decision-making in drug development process [49]. The toxicological analysis predicted that Pelargonidin and Cyanidin will be well-absorbed by the human intestine (HIA), which is a prerequisite for an orally administered potential drug candidate [49]. values as low as 0.77, 0.85 and 0.81 nM respectively. Thus, suggesting these compounds must be investigated in proper laboratory setup to validate these findings.

Conclusion:
The present study was designed to investigate the ERα binding property of Hibiscus sabdariffa polyphenols, largely the anthocyanis and to assess their toxicological risk. Accordingly, molecular docking study was performed to evaluate their potential to bind ERα active site (identical to the binding site of 17β-Estradiol).
All the anthocyanins formed comparatively stable conformations with ERα active site when compared to the control ligands, tamoxifen and raloxifen. These results can be further validated using a molecular dynamic simulation analysis. Further, with favorable toxicological parameters and potent bioactivity, all the anthocyanins exhibited strong affinity towards ERα, which supports the claim of estrogen-like activity of the plant extracts in MCF7 cell lines and in animal models. Among the anthocyanins, pelargonidin was the prominent find that satisfied all the in-silico analysis with promising outcome. These results will help in designing lead for novel ERα-agonists/antagonists as potential anticancer agents. Lastly, the ERα modulatory activity of these compounds must be validated in proper clinical or cell-based studies.

Declarations:
Funding: The authors did not receive support from any organization for the submitted work.

Class SMILE string Campesterol
Phytosterol    which provides the compound with good absorption.