Pharmacophore-Based Virtual Screening from Indonesian Herbal Database to Finding New Inhibitor of HDAC4 and HDAC7

Erlina, Andika, and Yanuar: Pharmacophore-Based Virtual Screening from Indonesian Herbal Database to Finding New Inhibitor of HDAC4 and HDAC7



Class IIa HDAC inhibitors have been developed as anti-cancer drugs, anti-inflammatory, Huntington’s disease, human papillomavirus and anti-diabetic. Compounds that have been developed into a class of HDAC inhibitors include hydroxamic acid, cyclic peptides, aliphatic acids and benzamide.1 However, these compounds are not specific because it can work as an HDAC inhibitor in class I and class II HDACs.2 Therefore, it is necessary to develop a class IIa HDAC inhibitor that is specific and potent to be used as a new drug candidate compounds.

The importance key of early level drug discovery is finding new lead compounds with potential interactions with the specific targets. Wet-lab high-throughput screening (HTS) is one of conventional method that adopted by pharmaceutical industry. High cost and low hit rate from HTS methods stimulate the researcher to develop computational methods (in silico) to conduct faster and low cost in drug discovery.3,4

Rapid search for small molecules for binding to its biological target become a key role in the drug discovery process. SBVS (structure-based virtual screening) is commonly used to facilitate virtual screening of a database, docking is one of virtual screening if 3D structure of target is available.5 Advances in structure biology by using X-ray crystallography and nuclear magnetic resonance spectroscopy can explain more detail about 3D structure of the targets and their interaction with ligands.6,7

The basic concept of a pharmacophore is an interaction between a ligand and a receptor as a function of an individual functional group’s chemical features. Since the chemical features in pharmacophore modelling are known and accepted by medicinal chemists, various models have been shown to be successful in the drug discovery process using a computational technique.8 Medicinal chemists use representations of pharmacophores to characterize the ability of molecular structures to bind to the site of specific biological targets. Since the ligand–receptor interaction takes place in the three-dimensional area, pharmacophore models that are based on observed 3D interaction patterns may be the better intuitional choice. Nevertheless, if a receptor-relevant ligand conformation or conformation ensemble is not yet known, the quantitative structure–activity relationship (QSAR) studies that are based on 3D pharmacophore models may diverge.9 One process of additional conformer generation, an often limiting and time-consuming step in pharmacophore matching methods, involves the 3D alignment of molecular features, for example, matching a screening molecule to a given pharmacophore model.10

Three dimension pharmacophore-ligand based use the database to find a compound hit more often used because it has several advantages, namely the complexity of the process of identifying a hit can be reduced by representation pharmacophoric of ligand interaction with a target that produces search time shorter and query-based pharmacophore can use to searching for a new drug candidate with a scaffold and different functional groups than the native ligand that is used for modeling pharmacophore.11

In recent studies, HDAC inhibitors have shown some essential pharmacophore features: (a) a functional group that acts as a chelating agent for the metal ion in the active site (anchor), (b) a hydrophobic cap for interaction with amino acid residues at the entrance of the N-acetyl lysine-binding channel, and (c) a linker with a 5–6 hydrocarbon chain length, which connects both the hydrophobic cap and the anchor.12,13 In this research, pharmacophore models were built from a dataset of known active ligands; the best-validated model would continue to the next step, virtual screening against herbald database contains 1377 compounds from Indonesian plant.14



Mac Mini with the specification Operating System X Yosemite version 10.10, processor 2.6 GHz Intel Core i5, memory 8 GB 1600 MHz DDR3, graphics Intel Iris 1536 MB.


The three dimensional (3D) sequence of protein HDAC Class IIa (HDAC4 and HDAC7) homo sapiens are retrieved from online database: NCBI website, Research Collaboratory for Structural Bioinformatics ( and website Protein Data Bank ( Two dimensional structures format .mol or .mol2 from database Indonesian Herbal database ( The three-dimensional structure of decoy are retrieved A Directory of Useful Decoys (DUD) (


Protein and dataset of ligand preparation

The three-dimensional structures of proteins HDAC4 and HDAC7 are obtained from PDB ( For HDAC4, 2VQM (PDB ID) structure of HDAC 4 catalytic domain bound to a hydroxamic acid inhibitor with resolution 1.8 Å.16 For HDAC7, 3C10 (PDB ID) crystal structure of catalytic domain of human histone deacetylase HDAC 7 in complex with Trichostatin A (TSA) with resolution 2.4 Å and 3ZNS (PDB ID) HDAC 7 bound with inhibitor Tfmo inhibitor TMP 942 with resolution 2.45 Å.17 A dataset of ligand is obtained from the MUBD-HDACs database.18 and divided into training and test set. Decoy compounds are generated from the test set using DUDE (

Formation and validation of 3D pharmacophore-ligand based models

Three-dimensional (3D) pharmacophore-ligand based are generated using LigandScout The dataset for HDAC4 includes 6 active compounds as a training set, 39 compounds as test set and 2040 decoys. The dataset for HDAC7 includes 11 active compounds as training set, 13 compounds as test set and 714 decoys. For each 3D pharmacophore-ligand based HDAC4 and HDAC7, 10 models of pharmacophore are build and copied into screening perspective. Test set and decoy are prepared into .ldb format file and used for screened against 10 pharmacophore models. Some parameters of validation such as ROC (Receiver Operating Characteristic), AUC100% (Area Under Curve), EF1% (Enrichment Factor) will be calculated.

Virtual screening

Based on validation of 3D pharmacophore-ligand based models, The best models were copied into screening perspective for virtual screening using LigandScout 4.09.1. In this step, herbaldb database ( that contains 1377 natural compounds original from Indonesian plants were screened with the best models from validation results and some new hits compounds would be obtained.14 Herbaldb database is prepared into .ldb format file.


Validation of 3D Pharmacophore-Ligand Based Models HDAC4

Based on validation parameter calculations in Table 1 and ROC graph (Figure 1), the model pharmacophore 3D-based ligands are best used in the model 6 and the model 10 with EF1% (24.0), the value AUC100% (0.65), accuracy (0.9793), sensitivity (0.3077), specificity (0.9922) and precision (42.86). On model of pharmacophore 6 and 10 there are seven features pharmacophore such as three HBA (Hydrogen Bond Acceptor) (red ball), one HBD (Hydrogen Bond Donor) (green ball), one aromatic ring (purple), one negatively ionizable area (red) and one hydrophobic (yellow ball) (Figure 2). The negatively ionizable area (-NH-) is a functional group that have interaction to the Zn2+ ion in the active site.

Figure 1

ROC best results of HDAC4 (pharmacophore model 6 and 10).
Figure 2

Pharmacophore features (2D and 3D) from the best model of HDAC4.
Table 1

Parameter validation of HDAC4 pharmacophore models.


Validation of 3D Pharmacophore-Ligand Based Models HDAC7

Based validation parameter calculations in Table 2 and ROC graph (Figure 3) model pharmacophore 3D-based ligands are best used in model 1 with a value EF1% (24.0), the value AUC100% (0.86), accuracy (0.5626), sensitivity (0.8462), specificity (0.5574) and precision (3.36). In the model, there are five features pharmacophore such as two HBA (red ball), one HBD (green ball), one negative ion (red) and one hydrophobic (yellow ball) (Figure 4). The negatively ionizable area (-NH-) is a functional group that have interaction to the Zn2+ ion in the active site.

Figure 3

ROC best results of HDAC7 (pharmacophore model 1).
Figure 4

Pharmacophore features (2D and 3D) from the best model of HDAC7.
Table 2

Parameter validation of HDAC7 pharmacophore models.


Virtual screening

Based on the election results and the best model for HDAC4 HDAC7 that has been done, then the next stage is pharmacophore models are used to perform virtual screening against a database of medicinal plants originally from Indonesia (herbaldb: totaling 1377 compound. From the results obtained HDAC4 pharmacophore screening models hit the compound number as many as 551 units (40.01% of the total in 1377 compounds). From the results obtained HDAC4 pharmacophore screening models hit the compound number as many as 647 units (46.99% of the total in 1377 compounds). Selection of hit compounds that will proceed to the stage of molecular dynamics simulations based on the examination of Drug-Likeness use of online applications molsoft ( Examination of drug-likeness includes several parameters Lipinski rule of five: a molecular weight that is < 500, the number of HBA is < 10, the number of HBD is < 5, and LogP values < 5.20 Compounds that are the hit artocarpesin, avicularin, dimboa glucoside, eriodictin, luteolin and mirabijalone C.

Based on Table 3, can be seen that the compounds hit DIMBOA glucoside value pharmacophore fit the highest of 44.60 with four pharmacophore features, three HBA and one HBD. Avicularin hit compounds, luteolin and artocarpesin have four features one hydrophobic pharmacophore, two HBA and one HDB. Compounds hit eriodictin, and mirabijalone c has three features of pharmacophore consisting of two HBA and one HBD.

Table 3

Hit compounds from virtual screening pharmacophore model against herbaldb.

CompoundsPharmacophore FeaturesPharmacophore-Fit Score
Mirabijalone C
Dimboa glucoside44.60


Based on Table 1 and 2, the best model for pharmacophore 3D-based ligands for HDAC4 is represented by model 6 and 10 which have score EF1% (24.0) and AUC100% (0.65) and the best model for pharmacophore 3D-based ligands for HDAC7 is represented by model 1 which have score EF1% (24.0) and AUC100% (0.86). The best Pharmacophore models are comply the requirements of enrichment factor indicator to define the enrichment performance for each matrix unit as very good (EFmax ≥ 30), good (30 > EFmax < 20), medium (20 > EFmax < 10), or poor (EFmax < 10).21 We analyzed based on pharmacophore model features of HDAC4 and HDAC7, there are four basic pharmacophore features that important for ligand-receptor interaction: HBA, HBD, negatively ionizable area and hydrophobic.

In the 2D pharmacophore of HDAC4, position HBA represented by O and OH groups, pharmacophore HBD group represented by N, pharmacophore negative ions represented by a group NH, hydrophobic and aromatic ring pharmacophore represented by benzene which is located to the right of the group O=S=O (Figure 2). In the 2D pharmacophore of HDAC7, position HBA represented by O and OH groups, pharmacophore HBD group represented by O, pharmacophore negative ions represented by a group NH, hydrophobic pharmacophore group represented by C in benzene marked in yellow (Figure 4).

From pharmacophore model of HDAC4 and HDAC7, CO-NH-OH known as Zinc Binding Group (ZBG) have two HBA, one HBD, and one negatively ionizable area; this group has an important role for interaction with the Zn2+ ion. Atom O from OH-NH-CO will build covalent coordination with Zn2+ ion, and the other atom will interact with active sites amino acids such as Asp and His.

Dimboa glucoside compound has one negatively ionizable area as Zinc Binding Group (ZBG) with the mechanism of interaction atom O from OH-NH-CO will build covalent coordination with the Zn2+ ion. The pharmacophore-fit score of dimboa glucoside (44.60) is the highest of all compounds. Artocarpesin, avicularin and luteolin have the pharmacophore-fit score (44.11), (44.40) and (44.33), they have one hydrophobic pharmacophore feature that similarly with the pharmacophore models features from HDAC4 and HDAC7. Mirabijalone c and eriodictin have the lowest pharmacophore fit-score (37.30) and (37.34) because they only have three basic pharmacophore features (2 HBA and 1 HBD).

All of the compounds have three basic pharmacophore features (2 HBA and 1 HBD), we are analyzed that addition of one negatively ionizable area (dimboa glucoside – Figure 5) will increase pharmacophore fit score rather than addition of one hydrophobic feature (artocarpesin – Figure 6, avicularin – Figure 7, and luteolin – Figure 8). Two compounds (mirabijalone c – Figure 9 and eriodictin – Figure 10) that only have three basic pharmacophore features (2 HBA and 1 HBD) have the lowest pharmacophore fit score. The mapping pharmacophore features from all compounds can be used further for chemical structure modification to add some missing pharmacophore features related to best pharmacophore models of HDAC4 and HDAC7.


Based on 3D pharmacophore-ligand based model validation for HDAC4 results, model 6 and 10 are the best models with seven pharmacophore features such as three HBA, one HBD, one aromatic ring, one negatively ionizable area and one hydrophobic. Based on 3D pharmacophore-ligand based model validation for HDAC7 results, model 1 is the best model with five pharmacophore features such as two HBA, one HBD, one negatively ionizable area and one hydrophobic. Based on virtual screening 3D pharmacophore-ligand based of HDAC4 and HDAC7 with herbaldb database results, there are six hit compounds (artocarpesin, avicularin, dimboa glucoside, eriodictin, luteolin and mirabijalone c). The results obtained in this research could be recommended for further studies such as molecular dynamic simulation, in vitro and in vivo.


Author (AY) thanks to the Publikasi Internasional Terindeks Untuk Tugas Akhir Mahasiswa (PITTA) 2016 grant from the Universitas Indonesia, for funding support.


The authors declare no conflict of interest.



Area Under Curve;


Directory of Useful Decoys;


Enrichment Factor;


Hydrogen Bond Acceptor;


Hydrogen Bond Donor;


Histone Deacetylase;


protein Data Bank;


quantitative structure –activity relationship;


Receiver Operating Characteristic;


Structure-Based Virtual Screening;


Trichostatin A;


Zinc Binding Group.



Cacabelos R , author. Epigenomic Networking in Drug Development: From Pathogenic Mechanisms to Pharmacogenomics. Drug Dev Res [Internet]. 2014;75(6):348–65.


Ferrari A, Fiorino E, Giudici M, Gilardi F, Galmozzi A, Mitro N , authors. et al. Linking epigenetics to lipid metabolism: focus on histone deacetylases. Mol Membr Biol [Internet]. 2012;29(7):257–66.


Ripphausen P, Nisius B, Peltason L, Bajorath J , authors. Quo Vadis, Virtual Screening? A Comprehensive Survey of Prospective Applications. J Med Chem [Internet]. 2010;53(24):8461–7.


Clark DE , author. What has virtual screening ever done for drug discovery? Expert Opin Drug Discov [Internet]. 2008;3(8):841–51.


Kroemer RT , author. Structure-based drug design: docking and scoring. Curr Protein Pept Sci. 2007;8(4):312–28


Villoutreix BO, Eudes R, Miteva M , authors. Structure-Based Virtual Ligand Screening: Recent Success Stories. Comb Chem High Throughput Screen. 2009;12(10):1000–16


Ghosh S, Nie A, An J, Huang Z , authors. Structure-Based Virtual Screening of Chemical Libraries for Drug Discovery. Curr Opin Chem Biol. 2006;10(3):194–202


Wolber G, Seidel T, Bendix F, Langer T , authors. Molecule-pharmacophore superpositioning and pattern matching in computational drug design. Drug Discov Today. 2008;13(1):23–9


Doweyko AM , author. 3D-QSAR illusions. J Comput Aided Mol Des [Internet]. 2004;18(7-9):587–96.


Yanuar A, Azminah EY, Andika EY, Erlina L, Syahdi RR , authors. In silico Approach to Finding New Active Compounds from Histone Deacetylase (HDAC) Family. Curr Pharm Des [Internet]. 2016;22(23):3488–97.


Schneider G , author. et al. “Scaffold-Hopping” by topological pharmacophore search: a contribution to virtual screening. Chem Int. 1999;38(19):2894–6


Kapustin GV, Fejér G, Gronlund JL, McCafferty DG, Seto E, Etzkorn FA , authors. Phosphorus-Based SAHA Analogues as Histone Deacetylase Inhibitors †. Org Lett [Internet]. 2003;5(17):3053–6.


Somoza JR, Skene RJ, Katz BN, Mol C, Ho JD, JC AJ , authors. Structural snapshots of human HDAC8 provide insights into the Class I histone deacetylases. Structure. 2004;12(7):1325–34


Yanuar A, Mun A, Bertha A, Lagho A, Syahdi RR, Rahmat M , authors. et al. Medicinal Plants Database and Three-Dimensional Structure of the Chemical Compounds from Medicinal Plants in Indonesia. Int J Comput Sci. 2011;8(5):180–3


Mysinger MM, Carchia M, Irwin JJ, Shoichet BK , authors. Directory of useful decoys, enhanced (DUD-E): Better ligands and decoys for better benchmarking. J Med Chem. 2012;55(14):6582–94


Bottomley MJ, Lo Surdo P, Di Giovine P, Cirillo A, Scarpelli R, Ferrigno F , authors. et al. Structural and Functional Analysis of the Human HDAC4 Catalytic Domain Reveals a Regulatory Structural Zinc-binding Domain. J Biol Chem [Internet]. 2008;283(39):26694–704.


Lobera M, Madauss KP, Pohlhaus DT, Wright QG, Trocha M, Schmidt DR , authors. et al. Selective class IIa histone deacetylase inhibition via a nonchelating zinc-binding group. Nat Chem Biol [Internet]. 2013;9(5):319–25.


Xia J, Tilahun EL, Reid TE, Zhang L, Wang XS , authors. Benchmarking methods and data sets for ligand enrichment assessment in virtual screening. Methods [Internet]. 2015;71:146–57.


Wolber G, Langer T , authors. Ligand Scout: 3-D Pharmacophores Derived from Protein-Bound Ligands and Their Use as Virtual Screening Filters. J Chem Inf Model [Internet]. 2005;45(1):160–9.


Lipinski CA , author. Lead- and drug-like compounds: the rule-of-five revolution. Drug Discov Today Technol [Internet]. 2004;1(4):337–41.


Huang N, Schoichet BK, Irwin JJ , authors. Benchmarking Sets for Molecular Docking. Journal of Medicinal Chemistry. 2006;49(23):6789–801. doi:10.1021/jm0608356. Benchmarking.