Bioinformatics: A rational combine approach used for the identification and in-vitro activity evaluation of potent β-Glucuronidase inhibitors

Identification of hotspot drug-receptor interactions through in-silico prediction methods (Pharmacophore mapping, virtual screening, 3DQSAR, etc), is considered as a key approach in drug designing and development process. In the current design study, advanced in-silico based computational techniques were used for the identification of lead-like molecules against the targeted receptor β-glucuronidase. The binding pattern of a potent inhibitor in the ligand-receptor X-ray co-crystallize complex was used to identify and extract the structure-base Pharmacophore features. Based on these observations; five structure-based pharmacophore models were derived to conduct the virtual screening of ICCBS in-house data-base. Top-ranked identified Hits (33 compounds) were selected to subject for in-vitro biological activity evaluation against β-glucuronidase enzyme; out of them, twenty compounds (61% of screened compounds) evaluated as actives, however eleven compounds were found to have significantly higher inhibitory activity, including compounds 1, 5–8, 10, 12–13, and 17–19 with IC50 values ranging from 1.2 μM to 34.9 μM. Out of the eleven potent inhibitors, seven compounds 1, 5, 6, 7, 8, 13, and 19 were found new, and evaluated first time for the β-glucuronidase inhibitory activity. Compounds 1, 5 and 19 exhibited a highly potent inhibition in uM of β-glucuronidase enzyme with non-cytotoxic behavior against the mouse fibroblast (3T3) cell line. Our combined in-silico and in-vitro results revealed that the binding pattern analysis of the eleven potent inhibitors, showed almost similar non-covalent interactions, as observed in case of our validated pharmacophore model. The obtained results thus demonstrated that the virtual screening minimizes false positives, and provide a template for the identification and development of new and more potent β-glucuronidase inhibitors with non-toxic effects.

Introduction β-Glucuronidase belongs to the glycosidase family of enzymes, which catalyze the hydrolysis of complex carbohydrates. The active site of the enzyme consists of a large cleft at the interface of two monomeric units. Two acidic amino acids, i.e., Glu 540 and Glu 451, and one aromatic amino acid, i.e., Tyr 504, have been proposed to be important for catalysis [1][2]. Human β-glucuronidase is homologous to the Escherichia coli enzyme β-glucuronidase. It catalyzes the hydrolysis of carbohydrates using two acidic a.a residues, Glu 540 and Glu 451. Additionally, the a.a residue Tyr 504 is involved in this catalytic event [3]. The catalytic mechanism involves three steps, as follows (1) Nucleophilic attack of the carboxylate anion on the anomeric carbon of sugar, (2) Hydrolysis of glucuronyl enzyme intermediate, and (3) De-glucuronidation [4][5].
Over-expression of β-glucuronidase enzyme activity is associated with several disorders, including various types of cancers, particularly hormone-dependent cancers, such as breast, prostate, and colon cancers. For the treatment of disorders associated with increased β-glucuronidase activity, D-saccharic acid 1, 4-lactone (DSL; saccharo lactone), silymarin, and silybin (crude drugs) are commercially available [6][7]. However, these drugs decreases immunity, and cause adverse effects. Therefore, there is a strong need to develop new β-glucuronidase inhibitors with improved potency and fewer adverse effects.
Structure-based pharmacophore mapping considered as a useful tool for medicinal chemists to identify novel ligands that have a high probability of being biologically active. This method utilizes the following steps: (I) Protein structure preparation, (II) Binding site detection, (III) Pharmacophore features identification, and (IV) Pharmacophore features selection.
Structure-based Pharmacophore can be efficiently used for virtual screening, ligand-receptor binding pose prediction, and binding site similarity search. Therefore, this method is a valuable tool for Hit and lead optimization, compounds library design, scaffold hopping, virtual screening, and multi-target drug design [8][9][10]. A successful virtual screening can identify molecules with novel chemical structural features that bind to the target receptor of interest in a large chemical space (e.g. the needle in a haystack concept).
The main purpose of our entire designed study was to develop structure-based Pharmacophore models, extracted with appropriate chemical structural features information, and use of those models to conduct the virtual screening of an in-house data-base in search of new lead candidates as inhibitors of β-glucuronidase with more potency [Fig 1]. For this purpose, we used advance in-silico techniques of computer-aided drug design (CADD) to reduce the large chemical space, and to increase the focus on more promising candidates towards lead discovery and optimization. model, 85 Hits were obtained by using 3LPG individual Pharmacophore model, and finally 100 Hits were obtained through the individual 3K4D Pharmacophore model; thus a total of 1,249 Pharmacophore-based Hits were identified [Fig 2].

Molecular docking results of Pharmacophore-based Hits
Molecular docking studies of pharmacophore-based Hits was carried out along with 66 previously reported inhibitors from literature (belonging to different classes of compounds) [13] [S1 Data-set] by using Fast Rigid Exhaustive Docking (FRED) software [14][15][16][17]. Rescoring of the chemgauss-4 scoring function of FRED software was performed by using GOLD docking software. Enrichment factors were calculated for 5%, 10%, 15%, and 20% of screened database for each scoring function chemgauss-4, chem score, gold score, and ASP score, respectively [ Fig 3A-3D] to examine the potential strength of all scoring functions for identifying drug-like candidates (redundancy of the in-house data-base), and to ultimately remove the non-binders (non-redundancy of the decoy set) [18], [S1 Appendix].
Receiver operating characteristic (ROC) curves. ROC curves are used to validate the docking software performance, which differentiates between the true binders (true positives) and non-binders (false positives). ROC curves are the plots between sensitivity and (1-specificity), sensitivity defines the presence of true actives in a data-base. Higher the sensitivity values Graphical representation of Pharmacophore-based Hits / non-Hits with decoys. The X-axis (blue bar-graph) represents total number of Pharmacophore-based Hits (1,249) compounds redundancy of in-house data-base, while the red bar-graph depicting total number of compounds non-redundancy in decoys. The Y-axis represents total number of screened compounds in the decoy-set (8,262) of in-house data-base.
https://doi.org/10.1371/journal.pone.0200502.g002 represent an increase number of true positives in the data-base, whereas specificity defines the presence of true-negatives (non-actives) in the data-base. The higher number of actives in the data-set represented the increased sensitivity and decreased probability of specificity (presence of non-binders). This statistical test is used in high-throughput computational-based virtual screening (HTS) for quick and efficient differentiation between the actives and non-actives in a decoy set of compounds. Based on this, we can consider an optimal model (true-binders) and reject (discard) sub-optimal model for data that is not required (comprised of nonbinders).
Accuracy can be measured by using area under the curve (AUC). An area of 1 represents a perfect test, whereas an area of 0.5 is supposed to be worthless; however, an area of 0.7-0.9 is considered as an acceptable value.
In our case study, the AUC value calculated as 0.76, which was a quite acceptable value, demonstrating that our data-set was considerably enriched with true binders (actives). Additionally, as we increased the cut off value of the data-base (small-subsets) in which true positives were present, the probability of finding TP increased, which ultimately raised the sensitivity [21][22], [Fig 4].

Structure-activity relationship (SAR)
SAR of halo (Cl, Br)-substituted indol derivatives. Halo-substituted bis-indols possess tremendous medicinal importance, [25] the docked pose interaction analysis of compound 8 demonstrated that H-bonding between the NH of indol and phenylalanine (Phe161) was responsible for its highly potent activity (44 folds higher than the standard D-saccharic acid 1, 4-lactone), along with π-π stacking interactions between the active site amino acid Tyr 472 and the phenyl ring of indol. The presence of di-hydroxy substituted phenyl in compound 8, one hydroxy group formed H-bonding with amino acid Lys 568. It was also responsible for the potent activity as compared to the activity of compound 17. This could be explained by the absence of one OH group, which was replaced with one OMe (methoxy) group, exhibited Harene interaction, a comparatively weaker interaction than the H-bonding. A decrease in the activity of compound 10 was observed due to the absence of one OH group, which was replaced with one OMe group, although this provided an electron donating group but it was not involved in H-bonding, as it was observed in compound 8. Glu 413 a.a also acted as a hydrogen bond acceptor for compound 9 within the binding pocket region of 5 Å. The decrease in the activity of compound 18 was also observed due to the absence of H-bonding of one OH substituent in the phenyl ring [Tables 2 and 3].
SAR of cyclic thioimidazole derivatives. Docked pose interaction analysis of phenyl sulfone-substituted cyclic thiourea derivatives provided an insight into the most active compound 7, from thioimidazole class which showed 78.1% inhibition,(IC 50 = 11.41 ± 0.04 μM). In this compound, NH formed H-bonding with the active site amino acid residues Glu 413 and Phe 161, whereas the di-methoxy phenyl substituent showed π-π stacking interactions with Tyr 472. On the other hand, compound 6 with 81.0% inhibition, and (IC 50 = 11.8 ± 0.86 μM), was slightly less potent than compound 7. This decrease in activity was likely due to the absence of H-bonding between NH and Phe 161 and π-π stacking interactions with Tyr 472 in compound 6. Compound 5 showed 95.1% inhibition, (IC 50 = 14.7 ± 1.88 μM) in comparison to the above two compounds. A lower activity was observed due to the absence of the di-methoxy substituent on the phenyl ring therefore the π-π stacking interactions did not exist. A decrease in the activity of compound 19 was also observed due to the absence of one OMe substituent on the phenyl ring; thus, the only phenyl ring possess less electron donating effect comparatively substituted with OMe, as a consequence the strength of the existing π-π stacking interactions becomes weaker, resulting to decrease in the biological activity [Tables 2 and 3].
https://doi.org/10.1371/journal.pone.0200502.g005 Amino acid Tyr 572 showed arene-arene, π-π stacking interactions with indol phenyl ring. Trp 549 acted as arene-H donor to the HB-acceptor methoxy oxygen for H-bonding. Thr556 behaved as HB-donor to the HB-acceptor lone pair of azo-nitrogen. NH of indol behaved as HB-donor to the HB-acceptor a.a Glu 413 for H-bonding within the binding region of receptor active site.

16.16±0.76
(Continued ) Discussion β-Glucuronidase is an important glycosidase enzyme with great biological importance, which provides space towards the libraries of small organic inhibitors designing [26][27]. It's over expression is associated with several cancers [28]. Therefore to overcome the adverse effects associated with the available drugs, there is a strong need to search and identify lead candidates which possesses the therapeutic potential against this target receptor with less adverse effects. For this purpose we first derived structure-based pharmacophore models using Ligand Scout software 3.0 version. These developed models were furthermore used to search and identify the drug-like candidates from in-house data-base by using MOE software. The Pharmacophore based features exist at similar distances in the respective screened compounds of data-base, to picked up (query editor) searched for those compounds which have possibly the similar Pharmacophore match with our derived models, (Pharmacophore search on the bases of distances b/w the respective features of functional-groups). For conducting the virtual screening, ICCBS in-house data-base was used as decoy set, which comprises over 9,000 structurally diversified molecules data-base was initially filtered to obtain drug-like candidates by using MOE software filter, removed all of those compounds which deviate from "Lipinski rule of five" and follow drug ability criteria, initially unwanted, highly reactive, toxic and those compounds which deviate from (ROF), and possess poor bio-availability, were removed from the data-base, the filters ultimately selected those compounds which possess good ADME / Tox properties. Software Omega (Open Eye Scientific Software) was used for energy minimization through its MM2 Molecular Mechanics force field parameters for each compound of data-base, atom types were corrected and molecular charges were also added through SYBYL software [29][30]. Virtual screening of filtered compounds identified small organic molecules, which follows drug-ability criteria described by Rule of five (ROF). Therefore, computational based highthroughput screening (HTS) helped out to remove or eliminate toxic, larger, unstable and non-drug-like candidates from data-base, in the way restrict to identify drug-like candidates. As a result it prevents to identify higher molecular weight compounds (usually transition metal complexes, peptides etc) which sometimes appear active at the in-vitro level but suffer with bio-availability issues at the in-vivo level [31][32].
The major purpose of in-silico based analyses is to examine the existing non-covalent interactions with catalytic amino acid residues through molecular docking binding poses which showed a good efficiency by using the software and good enrichment in the data-base with sufficient active compounds [33]. These knowledge-based research findings were furthermore validated with the help of statistical methods, including ROC curves, AUC values, and enrichment factors, which have been successfully applied to identify true positives (true-binders), and efficiently remove the false positives (non-binders) in the data-base. Therefore, structurebased Pharmacophore model can be efficiently used to quick search and identify the Pharmacophore features in various diverse classes of compounds and considered as an alternative tool of molecular docking, however this cannot be observed in the ligand-base Pharmacophore model which lacks the interaction with active site a.a residues, and usually Hits the similar type of compounds.
After in-silico based screening the top-ranked candidates were also examined using an invitro bio-assay screening for the evaluation of β-glucuronidase inhibitory activity [34]. Twenty compounds were selected, in which compounds 8 and 17 showed significantly higher inhibitory potential, whereas compounds 1, 5-7, 10, 12-13, 18-19, were also exhibited good  inhibitory effect however compounds 11, and 15 showed moderate inhibitory activity against β-glucuronidase enzyme, as compared to the standard, (D-saccharic acid 1, 4-lactone). Eleven potent compounds were further subjected to cytotoxicity assay against the mouse fibroblasts (3T3) cell line, to map their preliminary toxicity profile. The standard inhibitor used for 3T3 mouse fibroblasts cell line was cycloheximide. The potent leads of β-glucuronidase inhibitors 1, 5, and 19 exhibited completely non-cytotoxic behavior. Whereas compounds 6-8, 10, 12-13 and 17-18 showed a moderate cytotoxicity. The obtained results revealed that our potent inhibitors contained functional groups NH, sulfur, and thiourea. This observation was further supported from a recently published U.S. Patent [2012 / 0322797] December 20/2012 [35], which also included urea derivatives, sulfurcontaining compounds, methoxy-substituted quinoline derivatives, and halo-substituted phenyl-thiophene compounds. Moreover, we observed the presence of the phenyl ring, heteroatoms (sulfur and oxygen), NH (in case of urea), amide nitrogen, and oxygen atoms, which were all have important chemical structural features of the Pharmacophore, and have interacted with the amino acid residues of binding region, to establish important non-covalent interactions b/w ligand and receptor, including, H-bonding, H-arene, arene-arene, π-π stacking interactions, Van der waal forces of attraction and dipole-dipole interactions, these all contributed towards the best fit ligand-receptor binding pose and ultimately for potent biological activity.
In the way, bioinformatics computational based techniques have proved to be a useful rational approach towards drug designing and discovery process with improved potency and reduce toxic effects.

Structure-based individual pharmacophore mapping of PDB 3LPF
Two and three-dimensional (2D and 3D, respectively) structure-based Pharmacophore models were derived from PDB I.d 3LPF using Ligand Scout 3.0 version. The software was used to illustrate the following pharmacophore features (Fig 6A and 6B).
Two yellow spheres represented the hydrophobic methyl substituent; One yellow sphere showed one hydrophobic region due to aromatic substituent; One green vector (arrow head) depicted the H-bonding of NH with a conserved water molecule HOH 680; another green vector (arrow head) depicted the H-bonding of the hydroxyl group with Glu 413A.
Red vectors (arrow heads) showed the H-bonding of HOH 731 and HOH 733 with the hydroxyl group.

Structure-based individual pharmacophore mapping of PDB 3LPG
The 2D and 3D structure-based pharmacophore models were derived from PDB I.d 3LPG using Ligand Scout 3.0 version. The software depicted the following pharmacophore features (Fig 7A and 7B).
Two yellow spheres represented the hydrophobic region of methyl benzene ring interacting with the amino acids Val 473A, Mse 447A, and Phe 448A; Another yellow sphere showed hydrophobic aromatic ring and hydrophobic fluorine; the green vector (arrow head) depicted the H-bond donor NH for the acceptor HOH 667; The red vector (arrow head) represented the hydrogen acceptor of the carbonyl oxygen of aldehyde group from HB-donor HOH677. The 3D structure-based pharmacophore models were derived from PDB I.d 3K4D, highlighted with surface active site groove, depicted the features, two red vectors (arrow heads) for H-acceptors, one for the carboxylate anion, and one for the lactam carbonyl keto group. Three red-green vectors (arrow heads) showed the 3HB-donor / 3HB-acceptor of three OH substituent's group. One red pointed sphere depicted the negative ionizable area of carboxylate anion; along with the calculated distances in Å b/w the respective pharmacophore features. https://doi.org/10.1371/journal.pone.0200502.g008

Structure-based individual pharmacophore mapping of PDB 3K4D
The 2D and 3D structure-based pharmacophore models were derived from PDB I.d 3K4D using Ligand Scout software 3.0 versions. The software illustrated the following pharmacophore features (Fig 8A and 8B).
Two red vectors (arrow heads) represented HB-donor from a.a Tyr 472A, Arg 562A to HBacceptors, one is carboxylate anion (negative ionizable area), while another is lactam carbonyl keto group; three red-green vectors (arrow heads) depicted the 3HB-donor / 3HB-acceptor of three OH groups.

Generation of structure-based shared and merged feature pharmacophore models
A five-element structure-based shared feature pharmacophore model was also derived. For this purpose, we selected the PDB I.d 3LPF as a reference, to extract the keen binding pattern information's due to its potent inhibitor bound with IC 50 value of 326 nM. The remaining two individual derived Pharmacophore models were aligned over it in the alignment pane; finally a shared featured structure-based five element pharmacophore model was generated. Similarly, a merged feature model was also derived by extracting all the Pharmacophore features together. [Fig 9A and 9B].

Conclusions
In conclusion, our combined methodology of in-silico pharmacophore mapping and virtual screening established that computational virtual screening could provide a cost-effective and time saving approach for the selection of drug-like candidates. After experimental evaluation of the top ranked identified Hits, we obtained eleven potent inhibitors, three compounds with non-cytotoxic behavior while eight compounds with moderate cytotoxicity against the 3T3 mouse fibroblast cell line. In our designed case study virtual screening Hit results, (scaffold hopping) has been successfully performed, and we identified top 5% enriched data-base, contained new classes of compounds with potent biological inhibitory activity against the β-Glucuronidase enzyme.  9. (A) It depicts five element based shared feature pharmacophore model derived from PDB I.d, 3LPF, 3LPG and 3K4D by using Ligand Scout, three green vectors showed H-donors, five red vectors showed H-acceptors, four yellow spheres depicted lipophilicity with hydrophobic surface regions grey spheres showed the excluded volume, along with reference point set with a.a residues from active site contour.(B): It depicts a merged feature pharmacophore derived model, from PDB I.d 3LPF, 3LPG and 3K4D by using Ligand Scout, it comprised of the features, six red vectors showed the H-acceptors, four green vectors showed the H-donors, one red pointed sphere represented the negative ionizable area, six spheres showed the hydrophobic region, grey spheres depicted excluded volume along with reference point set with amino acid residues within active site contour.