Identification of a better Homo sapiens Class II HDAC inhibitor through binding energy calculations and descriptor analysis
© Tambunan and Wulandari licensee BioMed Central Ltd. 2010
Published: 15 October 2010
Human papillomaviruses (HPVs) are the most common on sexually transmitted viruses in the world. HPVs are responsible for a large spectrum of deseases, both benign and malignant. The certain types of HPV are involved in the development of cervical cancer. In attemps to find additional drugs in the treatment of cervical cancer, inhibitors of the histone deacetylases (HDAC) have received much attention due to their low cytotoxic profiles and the E6/E7 oncogene function of human papilomavirus can be completely by passed by HDAC inhibition. The histone deacetylase inhibitors can induce growth arrest, differentiation and apoptosis of cancer cells. HDAC class I and class II are considered the main targets for cancer. Therefore, the six HDACs class II was modeled and about two inhibitors (SAHA and TSA) were docked using AutoDock4.2, to each of the inhibitor in order to identify the pharmacological properties. Based on the results of docking, SAHA and TSA were able to bind with zinc ion in HDACs models as a drug target. SAHA was satisfied almost all the properties i.e., binding affinity, the Drug-Likeness value and Drug Score with 70% oral bioavailability and the carbonyl group of these compound fits well into the active site of the target where the zinc is present. Hence, SAHA could be developed as potential inhibitors of class II HDACs and valuable cervical cancer drug candidate.
The human papillomavirus (HPV) is a family of sexually transmitted, double-stranded DNA viruses with over 100 different genotypes identified till date. It is associated with many different types of cancers including cervical, vaginal, head and neck, penile and anal cancer. With approximately 450,000 newly diagnosed cases each year and a 50% mortality rate [1, 2], cervical cancer is the second most common cause of cancer-related death in women worldwide and it is almost always associated with HPV [3, 4]. Cervical cancer is the most common cancer of women in most developing countries, where it may account for as many as one fourth of female cancers .
HPV genotypes are divided into the low risk and high risk categories based on the spectrum of lesions they induce. The low-risk types induce only benign genital warts and include HPV 6 and 11. The high-risk group containing HPV 16, 18, 31, 33, 45 and 56 is associated with the development of anogenital cancers and can be detected in 99% of cervical cancers , with HPV16 found in 50% of cases .
A consistently effective and safe treatment for HPV infections is currently not yet available. Present therapeutic options are more directed at surgical eradication and/or by destroying malignant lesions via physical or chemotherapeutical intervention. A majority of these treatments have been developed empirically, few have been thoroughly tested, but none of them are completely satisfactory. In attemps to find additional drugs in the treatment of cervical cancer, inhibitors of the histone deacetylases have received much attention due to their low cytotoxic profiles .
The structural modification of histones is regulated mainly by acetylation/deacetylation of the N-terminal tail and is crucial in modulating gene expression, because it affects the interaction of DNA with transcription-regulatory non-nucleosomal protein complexes. The balance between the acetylated/deacetylated states of histones is mediated by two different sets of enzymes: histone acetyltransferases (HATs) and histone deacetylases (HDACs). HATs prefentially acetylate specific lysine substrates among other nonhistones protein substrates and transcription factors, affecting DNA-binding properties and, in turn, altering gene transcription. HDACs restore the positive charge on lysine residues by removing acetyl groups and thus are involved primarily in the repression of gene transcription by compacting chromatin structure. Therefore, open lysine residues attach firmly to the phosphate backbone of the DNA, preventing transcription. In this tight conformation, transcription factors, regulatory complexes, and RNA polymerases cannot bind to DNA Acetylation relaxes the DNA conformation, making it accessible to transcription machinery. High levels of acetylation of core histones are seen in chromatin-containing genes, which are highly transcribed genes; genes that are silent are associated with low levels of acetylation. Inappropriate silencing of critical genes can result in one or both hits of tumor suppressor gene inactivation in cancer .
Members of the classical HDAC family fall into two different phylogenetic classes, namely class I and class II. The class I HDACs (HDAC1, HDAC2, HDAC3, and HDAC8) are most closely related to the yeast (Saccharomyces cerevisiae). Class II HDACs (HDAC4, HDAC5, HDAC6, HDAC7, HDAC9, and HDAC10) share domains with similarity to HDA1, another deacetylase found in yeast .
The inhibition of HDAC activity by a specific inhibitor induces growth arrest, differentiation, and apoptosis of transformed cells as well as several cancer cells . Recent studies were directed to investigate the molecular effects of HDAC inhibition on cervical carcinoma cells as well as on primary human foreskin keratinocytes, separately immortalized with amphotropic retroviruses that carry the open reading frames of HPV 16 E6, E7 or E6/E7.
In these experiments one could show that E6/E7 oncogene function of human papillomavirus can be completely bypassed by HDAC inhibition. Both malignant and immortalized HPV 16/18-positive cells became blocked in G1/S transition despite ongoing viral gene expression. G1 arrest was accompanied by a down-regulation of cyclin D and cyclin A and a concomitant up-regulation of the cyclin kinase inhibitors (CKI) p21 and p27. Binding of both CKIs led to a complete block of the cyclin-dependent kinase (cdk2) activity and in turn prevented binding of E7. This was intriguing with respect to the reversibility of HPV transformation process, since it is thought that the abrogation of the growth inhibitory function of p21 and p27 through E7 represents a key event in HPV-induced carcinogenesis. HDAC inhibitors also trigger pRb degradation, while E2F expression remained unaffected. pRb degradation is an E7- specific phenomenon, since in E6-positive cells pRb only became hypophosphorylated. The presence of E2F under cell cycle arrest led to a classical "conflict situation" which finally induced apoptosis [11–13]. Hence, the knowledge how the transforming potential of HPV can be bypassed without switching off viral transcription could open new therapeutical perspectives for the treatment of cervical cancer .
The aim of this work is to analyze the interaction of Homo sapiens class II HDACs with SAHA and TSA that are already in the phase I/II clinical trials based on their binding affinity and pharmacological properties. Since, no theoretical works have been carried out in identifying the properties and specificity, we intend to identify the group that could act as potential binding inhibitors.
In this paper, we present homology models of six Homo sapiens Class II HDACs (HDAC4, HDAC5, HDAC6, HDAC7, HDAC9, and HDAC10) that are validated by comparison with the X-ray structure of HDAC4 and HDAC7, which became available during the course of our study. Two HDAC inhibitors (SAHA and TSA) are docked to the six homology models. The pharmacological properties of SAHA and TSA were identified using Molinspiration, Osiris Property, ToxBoxes, and Toxmatch-v1.06. Therefore, the molecular binding interactions between the histone deacetylases with SAHA and TSA were analyzed to provide some insights into the molecular interactions and designing new HDAC inhibitors.
Materials and methods
Sequence alignment and homology modeling
Sequences of HDAC4 (Genbank Accession Number NP_006028.2,410aa), HDAC5 (Genbank Accession Number NP_005465.2,383aa), HDAC6 (Genbank Accession Number NP_006035.2,264aa), HDAC7 (Genbank Accession Number NP_056216.2,422aa), HDAC9 (Genbank Accession Number NP_848510.1,376aa), and HDAC10 (Genbank Accession Number NP_114408.3,382aa) were extracted from the NCBI protein sequence database. The eleven three-dimensional structures of HDAC4 and HDAC7 (PDB code 2H8N, 2094, 2VQJ, 2VQM, 2VQO, 2VQQ, 2VQV, 2VQW, 3C0Y, 3C0Z, and 3C10) were used to find the catalytic domain of Homo sapiens Class II HDACs. All sequences were imported into the ClustalW program  and the sequence alignment editor, BioEdit, for multiple pairwise alignments. The resulting alignments were examined manually. The four dimensional structure of HDAC (PDB code: 1C3P, 1ZZ1, 2VQM, and 3EW8) was used as a template for Homo sapiens Class II HDACs homology modeling in the Modeller9v7 program . The multiple alignment between catalytic domain of class II HDACs and templates, using HHpred web server http://toolkit.tuebingen.mpg.de/hhpred to get an alignment pir format as an input format for homology modeling in the Modeller9v7 program .
Binding energy calculations
The Adaptive Poisson-Boltzmann Solver (APBS) software was used for binding energy calculations . Starting from the ligand conformations which is binding to Zn2+ ion in class II HDACs structure obtained by the docking process, were used as an input for binding energy calculations. The pdb format of class II HDACs and ligands complex was converted to pqr format by adding missing atoms, optimizing hydrogen bonding and assigning atomic charge and radius parameters. The CHARMM22 force field parameters were assigned to the class II HDACs. For ligand parameterization requires a MOL2-format representation of the ligand to provide the necessary bonding information. In this research, MOL2-format files were obtained through PRODRG web server http://davapc1.bioch.dundee.ac.uk/prodrg/.
Molecular descriptors calculation
Quantitative structure-activity relationships (QSARs) correlate the response with molecular properties of compounds under interest. Any compound to be considered as a lead must possess acceptable scores for all of the descriptors. Molinspiration http://www.molinspiration.com/ and Osiris Property explorer http://www.organic-chemistry.org/prog/peo/ were used to calculate twelve descriptors- logP, solubility, drug likeliness, polar surface area, molecular weight, number of atoms, number of rotatable bonds, volume, drug score and number of violations to Lipinski's rule for SAHA and TSA taken for the analysis .
Toxicity is defined as a dose-linked unsafe effect of a chemical compound on a target life form. Safety issues and ADME (Absorption, Distribution, Metabolism and Excretion) are major factors in drug failure and they are crucial to identify early in the discovery process to reduce late-stage attrition. ADME-Tox Box of PharmaAlgorithms http://pharma-algorithms.com/physchem.html was used to predict the various toxicity effects such as rat LD50, mouse LD50, oral bioavailability, Ames test, pKa and ion fraction values.
Results and discussion
Homology models validation
The final result of multiple sequence alignment of the class II HDACs is shown in Figure 2. HDAC4, HDAC5, HDAC6, HDAC7, HDAC9, and HDAC10 share the same catalytic pocket and zinc binding residues. Visualization of those residues using PyMOL, showed that area was solvent accessible and has negative charge electrostatic potential surface.
After sequence alignment of the Homo sapiens class II HDACs, five models for each of the six HDACs were generated through the Modeller9v7 program. The quality of the models thus obtained was then evaluated with respect to the conformation of the peptide backbone and the packing environment. Generation of the Ramachandran plots using VADAR (Volume, Area, Dihedral Angel Reporter)  gave no residues having Phi/Psi angles in the disallowed ranges, and the percentage of residues having Phi/Psi angles in the most favorable ranges is around 97-98%, similar to the template structures used. The quality of Ramachandran plots was satisfactory for all of the models. In the second test, the stereo/packing for residues of the same type in high quality experimental structures deposited in the Protein Data Bank were compared with our HDACs models using VADAR. A score of < 5.0 or worse usually indicate poor packing. In general, our HDACs models have similar packing scores compared to the template X-ray structure. A few residues are poorly packed in these models based upon scores less than 5, but those residues are not binding ligand domain. In summary, the quality of our HDACs models has been checked by two different criteria. The results showed that our models are reliable for performing further docking studies.
Comparison of the homology model and X-ray structure of HDAC4
During our homology modeling studies of HDACs, the X-ray structures of HDAC (PDB ID 1C3P, 1ZZ1, 2VQM, and 3EW8) cocrystallized with several inhibitors were reported at 1.80 Å, 1.57 Å, 1.80 Å, 1.80 Å resolution, respectively. Comparison of the X-ray structures to our homology model of Homo sapiens class II HDACs was an independent way of validating our result. With a backbone RMSD ≤ 0.5 Å for residues 648-1057 of HDAC4, 693-1075 of HDAC5, 579-842 of HDAC6, 521-942 of HDAC7, 633-1008 of HDAC9, and 40-421 of HDAC10, the Homo sapiens class II HDACs homology model is very similar to the X-ray structure 2VQM. The overall secondary and tertiary structures are very similar for the model and the X-ray structures.
Docking of SAHA and TSA to HDACs models
Free binding energy of SAHA and TSA - Homo sapiens Class II HDACs Complexes
Free Binding Energy
Homo sapiens Class II HDACs
Each ligand showed different affinities with the class II HDACs; for example, SAHA compound showed the best affinity with the HDAC 5 (-7.44 kcal/mol) based on AutoDock calculation and HDAC10 (-263.15 kJ/mol) based on APBS calculation. Whereas the same compound was found to be rank 2 with HDAC7 (-7.42 kcal/mol) based on AutoDock calculation and HDAC6 (-213.60 kJ/mol) based on APBS calculation, and rank 3 (-6.72 kcal/mol) with HDAC10 (AutoDock) and HDAC7 (-203.21 kJ/mol, APBS). Local free binding energy obtained from AutoDock of Homo sapiens class II HDACs complexed with an inhibitor showed that SAHA is a weaker inhibitor of HDACs than TSA. But, global binding energy of Homo sapiens class II HDACs and inhibitors obtained from APBS, showed that TSA to be a weaker inhibitor of HDACs than SAHA. There are differences in binding energy calculated with AutoDock and APBS; this is because AutoDock does not calculate columbic contribution from all of atoms in protein like APBS.
Molecular descriptors value for the SAHA and TSA
N Rotatable bonds
Toxicity risks of SAHA and TSA
Suberoyl anilide hydroxamic acid (SAHA)
Trichostatin A (TSA)
The ADME-TOX box results showed that the SAHA has an oral bioavailability of more than 70% i.e., good solubility and stability. It acts as a non-substrate and non-inhibitor of P-gp. SAHA does not undergo significant first-pass metabolism. Although the ability to inhibit P-gp is very important in cancer drug characteristics, but SAHA still has the possibility to become a candidate for cervical cancer drug, because cervical tissue is not a P-gp expressing tissue. The predicted LD50 values of SAHA in mouse and rat was found to be 1400 mg/kg and 1600 mg/kg respectively when administered orally. SAHA acts as a non-substrate when checked for the P-glycoprotein substrate specificity.
The three-dimensional models for six class II histone deacetylases (HDAC4, HDAC5, HDAC6, HDAC7, HDAC9, and HDAC10), which were built using homology modeling and validated by bioinformatics techniques and by comparison to an X-ray structures, were docked to SAHA and TSA. Our studies provide the structural view of the catalytic domain of a class II HDAC and reveal for this subclass specific features: (i) a novel zinc binding motif that is likely to be involved in substrate binding and/or protein-protein interactions and may provide a site for modulation of activity, and (ii) a unique active site topology in catalytic activity. SAHA and TSA predicted to inhibit the class II HDACs are effective to all forms HDACs. SAHA was satisfied almost all the properties i.e., binding affinity scores of SAHA in the six class II HDAC enzymes was -156.94 kJ/mol, -171.77 kJ/mol, -213.60 kJ/mol, -203.21 kJ/mol, -179.29 kJ/mol & -263.15 kJ/mol respectively, the Drug Likeness value (1.24) and drug score (0.37) with 70% oral bioavailability and the carbonyl group of these compounds fits well into the active site of the target where the zinc is present. Hence, SAHA could be developed as potential inhibitors of class II HDACs and valuable anti-cancer-agents.
This research is supported by University of Indonesia. We grateful to Prof. Teruna J. Siahaan, Ph.D., Department of Pharmaceutical Chemistry, The University of Kansas, Lawrence, USA for his excellent assistance in proofreading the manuscript and critical comments and to Dr. Ridla Bakri, Chairman of Department of Chemistry, Faculty of Mathematics and Natural Science, University of Indonesia for his support.
This article has been published as part of BMC Bioinformatics Volume 11 Supplement 7, 2010: Ninth International Conference on Bioinformatics (InCoB2010): Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/11?issue=S7.
- Parkin DM: The global health burden of infection-associated cancers in the year. Int J Cancer 2006, 118: 3030–3044. 10.1002/ijc.21731View ArticlePubMedGoogle Scholar
- Lowndes CM: Vaccines for cervical cancer. Epidemiol Infect 2006, 134: 1–12,. 10.1017/S0950268805005728PubMed CentralView ArticlePubMedGoogle Scholar
- Koutsky L: Epidemiology of genital human papillomavirus infection. Am J Med 1997, 102: 3–8. 10.1016/S0002-9343(97)00177-0View ArticlePubMedGoogle Scholar
- Kanodia S, Fahey LM, Kast WM: Mechanisms used by human papillomaviruses to escape the host immune response. Current Cancer Drug Targets 2007, 7: 79–89. 10.2174/156800907780006869View ArticlePubMedGoogle Scholar
- Kalvatchev Z, Rösl F: Human papilloma viruses: realities and perspectives. Biotechnol Biotechnol Equip 2007, 24: 137–143.View ArticleGoogle Scholar
- Walboomers JM, Jacobs MV, Manos MM, Bosch FX, Kummer JA, Shah KV, Snijders PJ, Peto J, Meijer CJ, Muñoz N: Human papillomavirus is a necessary cause of invasive cervical cancer worldwide. J Pathol 1999, 189: 12–19. 10.1002/(SICI)1096-9896(199909)189:1<12::AID-PATH431>3.0.CO;2-FView ArticlePubMedGoogle Scholar
- Bosch FX, Manos MM, Muñoz N, Sherman M, Jansen AM, Peto J, Schiffman MH, Moreno V, Kurman R, Shah KV: Prevalence of human papilomavirus in cervical cancer: a worldwide perspective. International biological study on cervical cancer (IBSCC) Study Group. J Natl Cancer Inst 1995, 87: 796–802. 10.1093/jnci/87.11.796View ArticlePubMedGoogle Scholar
- Thiagalingam S, Cheng KH, Lee HJ, Mineva N, Thiagalingam A, Ponte JF: Histone Deacetylases: Unique Players in Shaping the Epigenetic Histone Code. Ann N Y Acad Sci 2003, 983: 84–100. 10.1111/j.1749-6632.2003.tb05964.xView ArticlePubMedGoogle Scholar
- de Ruijter AJ, van Gennip AH, Caron HN, Kemp S, van Kuilenburg AB: Histone deacetylases (HDACs): characterization of the classical HDAC family. Biochem J 2003, 370: 737–49. 10.1042/BJ20021321PubMed CentralView ArticlePubMedGoogle Scholar
- Subha K, Kumar GR: Assessment for the identification of better HDAC inhibitor class through binding energy calculations and descriptor analysis. Bioinformation 2008, 3: 218–222.PubMed CentralView ArticlePubMedGoogle Scholar
- Finzer P, Kuntzen C, Soto U, zur Hausen H, Rösl F: Inhibitors of histone deacetylase arrest cell cycle and induce apoptosis in cervical carcinoma cells circumventing human papillomavirus oncogene expression. Oncogene 2001, 20: 4768–76. 10.1038/sj.onc.1204652View ArticlePubMedGoogle Scholar
- Finzer P, Ventz R, Kuntzen C, Seibert N, Soto U, Rösl F: Growth arrest of HPV-positive cells after histone deacetylase inhibition is independent of E6/E7 oncogene expression. Virology 2002, 304: 265–73. 10.1006/viro.2002.1667View ArticlePubMedGoogle Scholar
- Finzer P, Krueger A, Stöhr M, Brenner D, Soto U, Kuntzen C, Krammer PH, Rösl F: HDAC inhibitors trigger apoptosis in HPV-positive cells by inducing the E2F-p73 pathway. Oncogene 2004, 23: 4807–17. 10.1038/sj.onc.1207620View ArticlePubMedGoogle Scholar
- Acharya MR, Sparreboom A, Sausville EA, Conley BA, Doroshow JH, Venitz J, Figg WD: Interspecies differences in plasma protein binding of MS-275, a novel histone deacetylase inhibitor. Cancer Chemother Pharmacol 2006, 57: 275–81. 10.1007/s00280-005-0058-8View ArticlePubMedGoogle Scholar
- Thompson JD, Higgins DG, Gibson TJ: CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res 1994, 22: 4673–80. 10.1093/nar/22.22.4673PubMed CentralView ArticlePubMedGoogle Scholar
- Martí-Renom MA, Stuart AC, Fiser A, Sánchez R, Melo F, Sali A: Comparative protein structure modeling of genes and genomes. Annu Rev Biophys Biomol Struct 2000, 29: 291–325. 10.1146/annurev.biophys.29.1.291View ArticlePubMedGoogle Scholar
- Hildebrand A, Remmert M, Biegert A, Söding J: Fast and accurate automatic structure prediction with HHpred. Proteins 2009, 77(Suppl 9):128–32.View ArticlePubMedGoogle Scholar
- Dolinsky TJ, Czodrowski P, Li H, Nielsen JE, Jensen JH, Klebe G, Baker NA: PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res 2007, 35: W522–5. 10.1093/nar/gkm276PubMed CentralView ArticlePubMedGoogle Scholar
- Lipinski CA, Lombardo F, Dominy BW, Feeney PJ: Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev 2001, 46: 3–26. 10.1016/S0169-409X(00)00129-0View ArticlePubMedGoogle Scholar
- Willard L, Ranjan A, Zhang H, Monzavi H, Boyko RF, Sykes BD, Wishart DS: VADAR: a web server for quantitative evaluation of protein structure quality. Nucleic Acids Res 2003, 31: 3316–9. 10.1093/nar/gkg565PubMed CentralView ArticlePubMedGoogle Scholar
- Bottomley MJ, Lo Surdo P, Di Giovine P, Cirillo A, Scarpelli R, Ferrigno F, Jones P, Neddermann P, De Francesco R, Steinkühler C, Gallinari P, Carfí A: Structural and functional analysis of the human HDAC4 catalytic domain reveals a regulatory structural zinc-binding domain. J Biol Chem 2008, 283: 26694–704. 10.1074/jbc.M803514200PubMed CentralView ArticlePubMedGoogle Scholar
- Schuetz A, Min J, Allali-Hassani A, Schapira M, Shuen M, Loppnau P, Mazitschek R, Kwiatkowski NP, Lewis TA, Maglathin RL, McLean TH, Bochkarev A, Plotnikov AN, Vedadi M, Arrowsmith CH: Human HDAC7 harbors a class IIa histone deacetylase-specific zinc binding motif and cryptic deacetylase activity. J Biol Chem 2008, 283: 11355–63. 10.1074/jbc.M707362200PubMed CentralView ArticlePubMedGoogle Scholar
- Wang DF, Helquist P, Wiech NL, Wiest O: Toward selective histone deacetylase inhibitor design: homology modeling, docking studies, and molecular dynamics simulations of human class I histone deacetylases. J Med Chem 2005, 48: 6936–47. 10.1021/jm0505011View ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.