INTRODUCTION
Histone deacetylases (HDACs) are a class of enzymes that cleaves acetyl groups from e-N-acetylated lysine residues of histone proteins, which wrap DNA molecules. HDAC activity causes DNA molecules to be wrapped more tightly leading to various epigenetic regulations.1 In the last couple of decades, much effort was put into inhibiting HDAC activity as a strategy to design compounds against a wide range of conditions such as neurodegenerative diseases, inflammatory diseases, cancer, diabetes, cardiovascular diseases, and HIV infections.1,2,3,4
To date, 18 HDAC isoforms have been identified, which are classified into four classes (class I-IV) according to their sequence homology and tissue distribution (Table 1). Classes I, II, and IV are zinc-dependent classical HDACs, i.e., they require a Zn2+ ion in their catalytic site for activity, whereas class III, also known as sirtuins (SIRTs), are a structurally distinct class, which depends on NAD+ for catalytic activity.1 Currently, the crystallographic structures of the human HDAC1-4, 6-8, SIRT1-3, 5, and 6 are available, leading to an increase in high-throughput virtual screening (hVS) for design and identification of novel specific HDAC inhibitors5,6 as well as drug repurposing efforts.7 Drug repurposing has become a new approach in drug design as a means to reduce costs and attrition rates in clinical studies and speed up drug development process.8
Although isoform-selective HDAC inhibition is required in many HDAC-related treatment strategies9, pan-HDAC (e.g., trichostatin A, vorinostat, and valproic acid) and pan-SIRT inhibitors (e.g., nicotinamide) attract attention for their clinical effectiveness in diseases like cancer (Figure 1).10,11,12 Because most HDAC isoforms are associated with tumorigenesis and tumor progression, such polypharmacological approaches may prove more effective than isoform-specific inhibitors.13
Molecular docking is an in silico method to predict preferred binding orientation and affinity of a ligand with respect to a receptor. It has been used as a robust tool to identify hit matter as part of virtual screening for a long time. To improve virtual screening performance, consensus scoring is applied by combining scoring functions of multiple software programs, which usually is considered more accurate than single-score methods.6
In this study, we identified a number of drug molecules with potential to show pan-HDAC and pan-SIRT inhibitor activities using consensus structure-based hVS of the library of the Food and Drug Administration (FDA)-approved drugs in order to suggest drugs to be repurposed. The study also suggests pharmacophore hypotheses for virtual pan-HDAC and pan-SIRT drugs, through pharmacophore modeling, which are expected to improve our understanding of the design of potent de novo derivatives.
MATERIALS AND METHODS
Ligand preparation
The collection of FDA-approved drugs was obtained as 3D-coordinates sdf file from DrugBank (http://www.drugbank.ca) (accession: July 3, 2019).14 After removing the experimental, investigational, and nutraceutical compounds, the remaining 1502 ligands were prepared by removing salts and counter ions using LigPreP (2019-2, Schrödinger, LLC, New York, NY), and optimized geometrically using MacroModel (2019-2, Schrödinger, LLC, New York, NY) and conjugate gradients method according to OPLS_2005 forcefield parameters.15 The optimized library was directly used for Glide (2019-2, Schrödinger, LLC, New York, NY)16,17,18 and converted to sdf format for FRED (v3.3.1.2, Open Eye Scientific Software; Santa Fe, NM).19 For AutoDock Vina (v1.1.2, The Scripps Research Institute, San Diego, CA) the library was converted to pdbqt format by Open Babel (v2.4).20
Molecular docking protocol
The crystal structures of HDAC1-4, 6-8, and SIRT1-3, 5, 6 were downloaded from the RCSB Protein Data Bank (www.rcsb.org).21 The protein structures were prepared for docking by removing unwanted chains and residues, adding missing atoms, assigning hydrogen atoms, bond orders, partial charges (for Glide only), and setting ionization and tautomeric states, as well as the H bonds of the protein residues using the Protein Preparation Wizard of Maestro (Epik, Impact, Prime: 2019-2, Schrödinger, LLC, New York, NY).22 The prepared protein structures were assigned Gasteiger charges and converted to their pdbqt format using AutoDockTools (v1.5.7, The Scripps Research Institute, San Diego, CA) for AutoDock Vina. Grid maps of the active site of each protein were prepared using the receptor grid generation panel of Maestro (2019-2, Schrödinger, LLC, New York, NY) for Glide, Make Receptor (v3.3.1.2, Open Eye Scientific Software; Santa Fe, NM) for FRED. This procedure is automatically performed by AutoDock Vina. The central coordinates of the catalytic site of each structure were taken for a grid box of 27x103 Å3 size (see Electronic Supplementary Material for details). Molecular docking on Glide was performed at standard precision of 50 runs per ligand with the following settings: A scaling factor of 0.80 with 0.15 charge cut-off (absolute value) being applied for the ligands, Epik (2019-2, Schrödinger, LLC, New York, NY) state penalties were added to docking scores, nitrogen inversions and ring conformations were sampled, and post docking minimization was enabled. For FRED, docking was performed at high resolution mode with 50 runs per ligand. For AutoDock Vina, the default parameters were used and the virtual screening tool PyRx (v0.8, The Scripps Research Institute, San Diego, CA) was used to run the docking simulations on AutoDock Vina.23 Each ligand was assigned a docking score of the identified best pose from each software upon visual evaluation. The consensus score of a ligand for a given structure was determined by calculating the average of the three scores from the three software. A pan-HDAC and a pan-SIRT score were determined for each ligand by calculating the average of the consensus scores for all HDAC and SIRT structures, respectively.
Pharmacophore modeling
3D pharmacophore models for pan-HDAC and pan-SIRT inhibitor drugs were created with Phase (2019-2, Schrödinger, LLC, New York, NY)24 using the top-scoring three ligands according to each pan-HDAC (belinostat, bexarotene, and cianidanol) and pan-SIRT scores (alosetron, cinacalcet, and indacaterol) according to the following settings: Finding best alignment and common features method was applied, 50 conformers were generated for each ligand, three to five features were required in each hypothesis, all the query compounds were required to match each hypothesis, and the hypotheses (6 for pan-HDAC and 20 for pan-SIRT) were ranked according to PhaseHypoScore and the best hypothesis was selected for each group. The FDA-approved drug library was screened against each selected hypothesis using the Phase Ligand Screening panel with the default settings.
Statistical analysis was not performed in this study.
RESULTS AND DISCUSSION
Consensus structure-based VS
A total of 1502 FDA-approved drug molecules was in silico screened against HDAC1-4, HDAC6-8, SIRT1-3, 5, and 6 using three different docking software. For each drug molecule, a consensus score was assigned regarding each HDAC and SIRT isoform, which was the mean of the scores from the three software. A pan-HDAC score was then determined for each molecule by calculating the mean of the consensus scores for HDAC1-4 and HDAC6-8 (Table 2). The pan-SIRT scores were similarly calculated using the consensus scores for SIRT1-3, 5, and 6 of each drug molecule (Table 3). This approach assured that the molecules with good scores from all of the software and good consensus scores for all the isoforms ranked higher.
Virtual pan-HDAC drugs
Among the classical HDAC isoforms, classes I and II HDACs share a similar catalytic site topology. A zinc cofactor, chelating with two aspartate and one histidine residues, and a substrate, is at the bottom of a narrow lipophilic gorge that forms the catalytic site (Figure 2). Therefore, compounds with a linear lipophilic moiety with H-bond donor and acceptor groups at the tip such as trichostatin A can effectively occupy this site by chelating with the zinc and interacting with the zinc ligands. The active site of classical HDACs leaves little room for conformational flexibility, thus FRED, Glide, and AutoDock Vina usually produce similar poses for the molecules. Belinostat, bexarotene, and cianidanol were recorded as the top pan-HDAC scoring molecules in our study (Figure 3).
As expected for a pan-HDAC inhibitor anticancer drug, belinostat had the highest pan-HDAC score. Belinostat ranked 6th according to consensus HDAC1 scores and 1st according to consensus HDAC2, 6, and 8 scores (see Electronic Supplementary Material for details). Belinostat’s IC50 values of 7.20, 7.31, 7.82, and 7.16 nM have been reported against these targets.25,26,27,28 A hydroxamic acid at the end of a cinnamyl moiety is typical of potent classical HDAC inhibitors. The hydroxamic acid moiety of belinostat was in close contact with zinc and its ligands (Figure 4A-C). The cinnamyl benzene stacked with the histidine ligand (e.g., His709 of HDAC7) of zinc, as well as other aromatic sidechains of the nearby residues, e.g., Phe679 and Phe738 of HDAC7, in HDAC active sites (see Electronic Supplementary Material for details), corroborates findings of previous studies.29
Bexarotene, the second highest pan-HDAC scoring drug, is also an antineoplastic drug approved for the treatment of cutaneous T-cell lymphoma.30 The retinoid X receptor activator has not been tested against any HDACs so far, but it obtained the 7th best consensus HDAC2 and the 2nd best consensus HDAC6 score in our study. The benzoic acid moiety of bexarotene was mainly responsible for binding to HDAC active sites, in which the carboxylate group interacted with the zinc, its ligands, and nearby residues like Tyr308 of HDAC2. The benzene stacked with the aromatic side chains of the nearby residues such as Phe155 and Phe210 of HDAC2 (Figure 4D-F). Bexarotene was also the 8th best pan-SIRT scoring molecule, making it a likely inhibitor of the HDACs of all classes.
Cianidanol [(+)-catechin], a natural flavonol, has been withdrawn as a drug due to hematological toxicity. However, it is still marketed as an over-the-counter product for weight loss.31 It has been clinically evaluated for several cancer types but has no anti-HDAC activity record. In this study, it was the 2nd, 9th, and 4th best compound according to consensus HDAC4, 7, and 8 scores. Cianidanol is structurally different from the classical HDAC inhibitors regarding the zinc-interacting group, which is an ortho phenolic hydroxyl instead of a hydroxamic or carboxylic acid. Whereas these hydroxyl groups engaged with the zinc and its ligands, the benzene bearing the hydroxyls stacked with the aromatic side chains of the nearby residues such as Phe155 and Phe210 of HDAC2, as was recorded for belinostat and bexarotene (Figure 4G-I).
Virtual pan-SIRT drugs
The NAD+- dependent SIRT active site is composed of NAD+-binding Rossmann-fold subdomain and a distal zinc-binding pocket, and is more deeply buried and roomier than HDACs (Figure 2).32 The active pockets of SIRTs show large variations among the isoforms resulting in varied binding modes among the drug molecules and among the software for the same molecule. Since apo (NAD+ free) and holo (NAD+ including) states of SIRTs are both inhibited by inhibitors of diverse topology,32,33 we preferred the apo form of SIRTs in the hVS process to avoid missing out bulky drug molecules. Alosetron, cinacalcet, and indacaterol emerged as the best three drug molecules from the hVS study according to pan-SIRT scores (Figure 3).
Alosetron is a “setron”, i.e., 5-HT3 receptor antagonist, used for the treatment of irritable bowel syndrome.34 Although the effects of alosetron on SIRTs are yet to be studied, it obtained the best pan-SIRT score and the 7th best consensus SIRT2 score. Alosetron appeared to have two important moieties for interacting with the relevant key SIRT residues: (1) The imidazole ring which donates and accepts H bond (e.g., Arg71 of SIRT5) and makes π-π interactions (e.g., His158 of SIRT5), and (2) the indol-1-one that widely engages in π-π interactions (e.g., Tyr255 of SIRT5) (Figure 5A-C).
Cinacalcet is an allosteric activator of the calcium-sensing receptor and is used for the treatment of hyperthyroidism and hypercalcemia.35,36 The compound has no SIRT-related record but it obtained the 7th and 4th best consensus scores for SIRT1 and 2, respectively. The binding of Cinacalcet to SIRTs was supported by the π-π stacking via its two aromatic rings (e.g., His 58 and Tyr255 of SIRT5) as well as strong electrostatic contacts via the CF3 group with the active site residues (Figure 5D-F). For some SIRTs, the secondary amine formed H bonds (see Electronic Supplementary Material for details).
Indacaterol is a b adrenoceptor agonist used for the treatment of chronic obstructive pulmonary disease.37 The molecule, which ranked the third according to the pan-SIRT scores, has not been tested against SIRTs before. Indacaterol obtained the 2nd and 6th best consensus score for SIRT1 and 3, respectively. The two condensed ring systems of indacaterol greatly contributed to its theoretical binding affinity to SIRTs (Figure 5G-I). While these rings engaged in π-stacking, with residues like Phe64, His133, and Trp188 of SIRT6, the hydroxyl and amino groups made H bonds in most cases to further stabilize the binding (e.g., His133 and Leu186 of SIRT6).
Pharmacophore models for virtual pan-HDAC and pan-SIRT molecules
Pharmacophore models are hypothetical spatial orientations of the common pharmacophores (functional groups considered important for biological activity) for a set of ligands (or a single ligand) that share a biological property (activity, toxicity, etc.), which are widely exploited in rational drug design applications.38,39 We created a set of possible pharmacophore hypotheses for virtual pan-HDAC and pan-SIRT inhibitor drug molecules using top three scoring molecules of each group and selected the best hypothesis for each group according to PhaseHypoScore and BEDROC scores (scores showing how much a hypothesis fits to the query ligands in general) (Table 4).
The selected hypothesis for the virtual pan-HDAC inhibitor drugs (hypothesis 1) consists of a closely located H-bond acceptor (A) and donor (D) groups, and a distal ring (R). The alignment of belinostat, bexarotene, and cianidanol with hypothesis 1 suggests that A and D represent the hydroxamic, carboxylic, and phenolic groups interacting with the zinc and its ligands, while R aligns with the aromatic ring of these compounds that stack with the aromatic side chains of active site residues (Figure 6A-C). The best hypothesis for alosetron, cinacalcet, and indacaterol (hypothesis 2) comprises two adjacent rings for the condensed ring of these drugs, a vertical third ring regarding the two for a separate aromatic group, and a hydrophobic group (H) close to the third ring representing a hydrophobic substitution to the third ring, namely methyl, trifluoromethyl, and ethyl (Figure 6D-F). Thus, hypothesis 2 reflects the hydrophobic nature of SIRT catalytic site. Cianidanol and cinacalcet showed the best alignment to hypothesis 1 and 2, respectively (see Fitness score in Table 4).
To test the accordance of these hypotheses with the consensus structure-based hVS campaign, we screened the drug molecules against both hypotheses and compared the results with respect to pan-HDAC and pan-SIRT scores from the molecular docking by calculating the 10% enrichment factor (Table 4). This metric shows how many of the drug molecules that are among the top 150 pan-HDAC scoring molecules (i.e., top 10% drugs) are listed in the top 150 compounds according to PhaseScreenScore (a score that shows how much a screened molecule fits to the pharmacophore hypothesis) for hypothesis 1. The same applies for pan-HDAC scores and hypothesis 2. The 10% enrichment factor was 20% for both hypotheses, showing that both methods predicted somewhat similar drug molecules as top virtual pan-HDAC and pan-SIRT inhibitors.
CONCLUSION
A total of 1502 FDA-approved drugs were screened against a set of classical HDACs and SIRTs with available crystal structures using FRED, Glide, and AutoDock Vina. The drug molecules were ranked according to their average consensus HDAC and SIRT scores to identify the drug molecules that can potentially inhibit as many HDAC or SIRT isoforms, i.e., virtual pan-HDAC and pan-SIRT inhibitors. The consensus approach in this method works in two ways: Consensus among the software used and among the isoforms. Belinostat, bexarotene, and cianidanol were the best scoring virtual pan-HDAC inhibitors. Although bexarotene and cianidanol have not been tested against HDACs, they have potential against HDACs and could be repurposed for HDAC-related indications. Specifically, bexarotene may show potent in vitro activity against HDAC2 and 6; and cianidanol against HDAC4, 7, and 8. Among these molecules, belinostat is already a confirmed pan-HDAC inhibitor used as an anticancer agent, which shows the effectiveness of the hVS methodology. On the other hand, other known pan-HDAC inhibitors such as vorinostat were not listed among the top-scoring drugs suggesting that the hVS method has weaknesses as well. Alosetron, cinacalcet, and indacaterol obtained the best pan-SIRT scores. Although these compounds have no SIRT record, they may be useful in pan-SIRT-related conditions. Alosetron could be a promising inhibitor of SIRT2, cinacalcet of SIRT1 and 2, and indacaterol of SIRT1 and 3. Bexarotene was also listed among the top-ten pan-SIRT scoring drugs, which could be a potent inhibitor of all HDAC classes and repurposed for a unique indication in this regard. Taken together, these drug molecules may find new indications related to pan-HDAC and pan-SIRT inhibition.
Two pharmacophore hypotheses were formulated, one for top pan-HDAC scoring drug molecules (hypothesis 1) and the other for pan-SIRT scoring drug molecules (hypothesis 2). The top three pan-HDAC and pan-SIRT scoring drugs aligned very well with their respective hypothesis. The pharmacophore features in these hypotheses were in compliance with the binding interactions of the drug molecules predicted by docking. Rankings of the drug molecules according to the pharmacophore hypotheses and molecular docking screens bore similarities, i.e., some of the top-scoring molecules in pharmacophore screens were also among the top pan-HDAC and pan-SIRT scoring drugs. Therefore, the ligand- and structure-based hVS methods yielded compatible results.
Click for Electronic Supplementary Material for details access link: https://cms.galenos.com.tr/SolvePark/Uploads/Files/1e021a9828a942e8a2176106560fb2df.pdf
Conflict of interest: No conflict of interest was declared by the authors. The authors are solely responsible for the content and writing of this paper.
Funding sources:This study was funded by a grant from Scientific and Technological Research Council of Turkey (TÜBİTAK) (grant no: 115S387).