Selection of an improved HDAC8 inhibitor through structure-based drug design.

Histone deacetylases (HDACs) are enzymes, which catalyze the removal of acetyl moiety from acetyl-lysine within the histone proteins and promote gene repression and silencing resulting in several types of cancer. HDACs are important therapeutic targets for the treatment of cancer and related diseases. Hydroxamic acid inhibitors show promising results in clinical trials against carcinogenesis. 120 hydroxamic acid derivatives were designed as inhibitors based on hydrophobic pocket and the Zn (II) catalytic site of HDAC8 active site using Structure Based Drug Design (SBDD) approach. High Throughput Virtual screening (HTVs) was used to filter the effective inhibitors. Induced Fit Docking (IFD) studies were carried out for the screening of eight inhibitors using Glide software. Hydrogen bond, hydrophobic interactions and octahedral coordination geometry with Zn (II) were observed in the IFD complexes. Prime MM-GBSA calculation was carried out for the binding free energy, to observe the stability of docked complexes. The Lipinski's rule of five was analyzed for ADME/Tox drug likeliness using Qikprop simulation. These inhibitors have good inhibitory properties as they have favorable docking score, energy, emodel, hydrogen bond and hydrophobic interactions, binding free energy and ADME/Tox. However, one compound (Cmp22) successively satisfied all the studies among the eight compounds screened and seems to be a promising potent inhibitor against HDAC8.

Background

Histones are small basic proteins with two domains. The Cterminal domain is located inside the nucleosome core and the N-terminal domain is rich in lysine residues extending out of the nucleosome core [1]. Positive charge of the lysine residues makes the histones to bind with the negatively charged DNA that mainly regulates the interactions between DNA and histones. A single nucleosomal unit consists of octameric histones (H3–H4 tetramer and H2A–H2B dimer) coiled by a stretch of 146 base pairs of DNA. Among the above four histones, H3 and H4 are mainly targeted for various posttranslational modifications including acetylation, phosphorylation, and methylation [1, 2]. Among such modifications, in particular, acetylation has been highlighted as an important mechanism in post-transcriptional determination [2, 3]. Acetylation of the histones can be controlled by two counteracting enzymes in the form of antagonist such as Histone acetyl transferases (HATs) and HDACs. Addition and removal of acetyl moieties were catalyzed from the ε -amino groups of lysines [4-6]. HDACs are important gene expression regulators in eukaryotic cells affecting angiogenesis, cell-cycle arrest, apoptosis, terminal differentiation of different cell types [7, 8]. Phylogenic analyses reveal human HDACs into four distinct classes, Class I-IV. Class I includes HDACs 1-3 and 8 that are homologous to yeast Rpd3 (exclusively expressed in the nucleus). Class II includes HDACs 4-7, 9 and 10 that are homologous to yeast Hda1 (shuttled between the cytoplasm and the nucleus). Class III includes Sirtuins 1-7 that are homologous to yeast sir2. Class IV, the recently identified one includes HDAC 11. Class III HDACs are structurally unrelated to the other three classes of HDACs and Class IV HDACs possess structural homology and has only a low overall sequence similarity with both classes I and II. HDACs 1-11 are metalloenzymes that require Zn (II) at catalytic site for deacetylation mechanism. HDACs in classes I and IV are 42-45 kDa, whereas class II HDACs are about 120-130 kDa in size [9-11]. The malignancy of leukemias and lymphomas are linked to the aberrant recruitment of HDACs [12, 13]. The HDAC inhibitors (HDACIs) are reported to have anti-tumor activity in preclinical models and in clinical trials, promisingly to become new antineoplastic drugs [14]. HDACIs are shown to alter gene transcription and exert various anti-tumor effects [11, 15, 16]. HDAC8 protects the human ever-shorter telomerase 1b protein, a telomerase activator from ubiquitin-mediated degradation [13], associates with R-actin cytoskeleton and plays a vital role in the regulation of contractility in differentiating smooth muscle cells [14]. Recent studies indicate that HDAC8 inhibitor induces apoptosis in T cell derived tumor cells but not increasing the histone or tubulin acetylation levels [11]. Suberoylanilide hydroxamic acid (SAHA), a HDAC inhibitor has been clinically validated and approved by Food and Drugs Administration in the treatment of cutaneous T cell lymphoma as a new therapeutic drug [17]. In the recent computational studies [18-21], it was observed that the known HDAC inhibitors and newly designed HDAC8 inhibitors docked at the active site of HDAC8 are showing the potent results. These results could be used to synthesize the HDAC8 inhibitors or suggesting designing some novel inhibitors against HDAC8 in further studies. Our current study aims to use the SBDD approach for designing new compounds which are structurally related to SAHA inhibitor. Simultaneously, HTVs, IFD and Prime MM-GBSA studies were also carried out. Then, hydrogen and hydrophobic interaction studies were analyzed in the docked complexes. The results of in silico and interactions studies of SBDD complexes were compared with that of in silico results of SAHA to show the inhibitory potency of the designed inhibitors. These comparisons suggest that the designed SAHA like compounds can be used as HDACIs against carcinogenesis in future. The designed HDACIs may be potent in the in vitro and in vivo studies.

Methodology

All the molecular modeling calculations like rigid (HTVs) and flexible (IFD) docking, binding free energy calculation, and ADME/Tox were carried out using Glide software (Schrödinger LLC 2009, USA) in CentOS EL-5 workstation. PyMOL software was used for graphical visualization, analysis of the hydrogen bond interactions and Zn (II) coordination geometry to produce quality images. Hydrophobic interactions were observed between the active site of target and ligand using Ligplot software [22].

Protein Preparation

The X-ray crystal structure of human HDAC8 complexed with SAHA (1T69) was recovered from Protein Data Bank (PDB: www.rcsb.org). The deposited structure of PDB was modified viz addition of hydrogen atoms, assigning correct bond orders, fixing of the charges and orientation of groups for Glide modeling calculations. The orientation of hydroxyl groups and amide groups of ASN, GLN and HIS were converted into the charged state. The amino acid flips were assigned and H-bonds optimized iteratively. Non-hydrogen atoms were energy minimized until the average root mean square deviation reaches 0.3Å. Schrödinger modules Glide, Prime, QSite, Liasion and MacroModel were used for protein preparation.

Structure Based Drug Design and Ligand Preparation

SBDD has the potential significance in pharmaceutical related researches, especially in the new drug design progress [23-24]. Application of SBDD techniques is supported by an exponential growth in the number of experimental protein 3D structures from X-ray crystallography or NMR spectroscopy [23, 25]. The HDAC8 active site consists of a long, narrow and primarily hydrophobic tunnel formed by G151, F152, H180, F208, M274, and Y306. The end of the hydrophobic tunnel contains a catalytic site [26]. Designing a series of 120 hydroxamic acid derivatives (data not shown) based on the active site residues and the catalytic Zn (II) using SBDD approach has been carried out for the molecular modeling studies. The basic core of the molecules is similar to SAHA inhibitor. Designed inhibitors were subjected into ligand preparation by Ligprep 2.3 module (Schrödinger). Ligprep performs addition of hydrogens, 2D to 3D conversion, realistic bond lengths and bond angles, low energy structure with correct chiralities, ionization states, tautomers, stereochemistries and ring conformation.

High Throughput Virtual screening

HTVs is one of the methods to screen one or more compounds from a set of compounds. Glide performs HTVs that requires previously calculated receptor grid and one or more ligand structures. The active site of the receptor and its properties were calculated on a grid that provides accurate scoring function along with energy when ligand is docked. The receptor grid allows the ligands to bind in more than one possible and meaningful conformation. During the receptor grid generation process, the minimized crystal structure bound with SAHA was loaded in the workspace; the active site of the receptor was calculated automatically and shown as a grid box by picking the ligand. In HTVs, the flexibility of the receptor is restricted, but van der Waals radii of non-polar atoms were calculated to decrease close contact penalties between ligand and active site residues. All the designed compounds were screened against the calculated grid box using Standard Precision (SP) docking and 30 best compounds were selected based on the Docking score, Glide energy, Glide emodel and non-bonded interactions. The selected compounds were subjected to Extra Precision (XP) docking, which is a more powerful scoring and discriminating procedure, where the receptor was held rigid while docking. XP docking took more time for screening compared to SP docking. Eight compounds (Figure 1) were chosen by the above selection method to proceed with IFD.

Figure 1

Chemical diagrams of hydroxamic acid derivatives used in the study.

N1-hydroxy-N8-phenyloctanediamide [SAHA],

N1-(7-(9H-carbazol-4-yl)-1H-indol-3-yl)-N8- hydroxyoctanediamide [Cmp14],

N1-(4-(9H-carbazol-1-yl)- 1H-indol-3-yl)-N8-hydroxyoctanediamide [Cmp17],

N1-(6-(9H-carbazol-3-yl)-1H-indol-3-yl)-N8-hydroxyoctanediamide [Cmp22],

N1-(6-(2-amino-2-formylethyl)-1H-indol-3-yl)-N8-hydroxyoctanediamide [Cmp36],

N1-hydroxy-N8-(5-phenyl-1H-indol-3-yl) octanediamide [Cmp52],

3-(N1-hydroxyoctanediamido)-1H-indole-6-carbonyl chloride [Cmp108],

N1-(5-formyl-1H-indol-3-yl)-N8-hydroxyoctanediamide [Cmp118],

N1-(6-formyl-1H-indol-3-yl)-N8-hydroxyoctanediamide [Cmp119].

Induced Fit Docking

In SP and XP docking studies, ligands were docked at the active site of a receptor, where the receptor was held rigid and the ligands were free to rotate. However, rigid receptor with flexible ligand may provide misleading results in docking simulation. When ligand binds at the active site, it undergoes side-chain or backbone conformation or both in many proteins. These conformational changes allow the receptor to generate close conforms to the shape and binding mode of the ligand. This is known as IFD. IFD is one of the main complicating factors in SBDD studies and predicts accurate ligand binding modes and concomitant structural movements in the receptor using modules of Glide and Prime. Python scripts automate the IFD process and an interface in Maestro (Schrödinger). Energy minimized human HDAC8 protein complexed with SAHA was loaded in the workspace and the ligand was selected to specify the active site. Then, the IFD calculations were carried out for eight compounds (from SP and XP screening). Van der Waals radii of non-polar atoms of the receptor and ligand were scaled by a factor of 0.50, and 20 conformational poses were calculated. Based on the docking score, Glide energy, Glide emodel and non bonded interactions noticed in IFD results, the best conformation was chosen for further calculations.

Prime MM-GBSA

Prime MM-GBSA approach was used to calculate ligand binding energies and ligand strain energies for a ligand and a single receptor. MM-GBSA is a method that combines OPLSAA molecular mechanics energies (EMM), an SGB solvation model for polar solvation (GSGB), and a non-polar solvation term (GNP) composed of the non-polar solvent accessible surface area and van der Waals interactions. Here, the IFD pose viewer file of the best conformation chosen was given as the source in Prime MM-GBSA simulation. The total free energy of binding: ΔGbind = Gcomplex – ( Gprotein + Gligand ), where G = EMM + GSGB + GNP

ADME/Tox

QikProp, the prediction program designed by Prof. William L. Jorgensen was used to calculate ADME (Absorption, Distribution, Metabolism, and Excretion) properties. Qikprop is quick, accurate and predicts physically significant descriptors and pharmaceutically relevant properties of organic molecules. Ligprep minimized ligands were given as a source in Qikprop 3.2. Qikprop modules provide the ranges of molecular predicting properties for comparing the properties of a particular molecule with those of 95% of known drugs (Schrödinger 2009).

Discussion

SAHA was redocked at the active site of the target. The docked complex has a docking score of ­2.37, glide energy ­44.79 and glide emodel ­50.16 (Table 1 see supplementary material). The ligand maintains four hydrogen bond interactions with the active site residues, namely, D178 (OD2), H143 (ND1), H180 (NE2) and Y306 (OH) (Table 2 see supplementary material). The Zn (II) ion holds six coordination interactions with the catalytic site residues. In the octahedral geometry formed, it maintains four interactions with the active site residues and two interactions with O1 and O2 atoms of SAHA. The coordination distance range is 2.04­2.47 Å (Table 3 see supplementary material). SAHA has hydrophobic interactions with G151, F207, and F208 residues at the active site (Table 4 see supplementary material). Compound (Cmp) 14 docked complex orientation possesses good docking score ­6.25, glide energy ­62.69 and glide emodel ­87.81 (Table 1). Cmp14 maintained three hydrogen bond interactions with D178 (OD2), H180 (NE2) and Y306 (OH). In addition, two hydrogen bond interactions were observed with G140 (O) and Q263 (NE2) residues (Table 2). Zn (II) ion forms octahedral coordination with certain amino acid residues and the ligand. Five coordination interactions were formed with the active site residues and one interaction with O2 of the Cmp14. The coordination distance range is 2.01Å­2.37Å (Table 3). Six hydrophobic interactions were noticed in the complex viz. H143, F152, C153, F207, F208 and M274 (Table 4).

Cmp 17 docked at the active site of target has good docking score ­6.64, glide energy ­55.09 and glide emodel ­76.13 (Table 1). Interestingly, the ligand has seven hydrogen bond interactions with active site residues H143 (NE2), D178 (OD2), G206 (O), F208 (N), D267 (OD2), G303 (O) and Y306 (OH) (Table 2). Zn (II) ion forms octahedral coordination geometry with the catalytic site residues and the ligand. The ligand atoms O2 and O3 retained two coordinations with Zn (II) ion and active site residues have four, all in the distance range 2.07­ 2.49Å (Table 3). Cmp17 maintains four hydrophobic interactions with G151, F152, H180, and M274 residues at the tunnel-shaped hydrophobic region (Table 4). For Cmp 22, docking score is ­7.00, glide energy is ­67.02 and glide emodel is ­95.72 (Table 1). The ligand exhibits five hydrogen bond interactions, three interactions with H180 (NE2), G206 (O) and G303 (O) and remaining two hydrogen bonds with Y306 (OH) (Table 2 & Figure 2a). Zn (II) ion is in octahedral coordination, five with active site residues and remaining one with the ligand. Zn (II) coordination distances range is in the range 2.00­2.38Å (Table 3 & Figure 2b). H143, G151, F152, F207, and F208 have hydrophobic interactions with the ligand at the active site (Table 4). Cmp 36 docked complex has a docking score of ­6.38, glide energy of ­61.09 and glide emodel of ­82.68 (Table 1). The ligand has nine hydrogen bond interactions with G140 (O), D178 (OD2), H180 (NE2), K202 (NZ), G206 (O), Q263 (NE2), M274 (O), S276 (OG) and Y306 (OH) residues (Table 2). Zn (II) ion forms five coordinations with active site residues and one with the ligand. Octahedral coordination is observed and distances range is 2.01­2.30 Å (Table 3). The ligand is involved in five hydrophobic interactions with H143, G151, F152, F207 and F208 (Table 4). For Cmp 52, five hydrogen bond interactions with active site residues Y100 (OH), G140 (O), D178 (OD2), G206 (O) and Q263 (NE2) (Table 2) are observed. Ligand has favorable docking score of ­5.93, glide energy of ­56.54, and glide emodel of ­76.55 (Table 1). Zn (II) ion maintains octahedral coordination geometry, one with the ligand and remaining five with active site residues. Coordination distance range is 2.06­2.42 Å (Table 3). This ligand has hydrophobic interactions with H143, G151, F152, F208 and M274 residues (Table 4).

Figure 2

a) Hydrogen bond interactions of Cmp22 with active site residues; b) Octahedral coordination geometry of Cmp22 with the catalytic site residues.

Cmp 108 docked complex has a docking score ­6.25, glide energy ­53.99, and glide emodel ­83.37 (Table 1). The ligand has six hydrogen bond interactions with the active site residues G140 (O) (two donors), D178 (OD2), H180 (NE2), Q263 (NE2) and Y306 (OH) (Table 2). The ligand has only one co-ordination and the active site residues have five coordinations with Zn (II) ion. The octahedral coordination distances range is 2.01­2.30 Å (Table 3). H143, G151, H180, F273, M274 residues have hydrophobic interactions with the ligand (Table 4). Cmp 118 docked at the active site shows a docking score ­6.25, glide energy ­62.69, and glide emodel ­87.81 (Table 1). Four hydrogen bond interactions of the ligand were noted with the active site residues K33 (O), G140 (O), D178 (OD2), and Y306 (OH) (Table 2). Zn (II) ion forms octahedral coordination with active site residues and ligand. The distance range is 2.00­3.25 Å (Table 3). Seven hydrophobic interactions are found with H143, G151, F152, C153, F208, M274 and Y306. Y306 has a hydrogen bond interaction with O3 of the ligand, coordination formed with Zn (II) ion and hydrophobic interactions with some of atoms in ligand (Table 4).

Cmp 119-bound complex has a docking score ­7.33, glide energy ­57.90 and glide emodel ­87.49 (Table 1). The active site residues H142 (NE2), D178 (OD2), H180 (NE2) (two donors), K202 (NZ), G206 (O), D267 (OD2) and Y306 (OH) have a total of eight hydrogen bond interactions with ligand (Table 2). In the octahedral coordination of the amino acids and ligand atoms with Zn (II) the distance range is 1.99­3.23 Å (Table 3). H143, G151, F152, F207, F208 and M274 residues show six hydrophobic interactions with the ligand (Table 4).

Prime MM-GBSA and ADME/Tox

Prime MM-GBSA calculation was carried out for SAHA as well as for the eight compounds designed by SBDD approach. SAHA has the binding free energy (ΔG) of ­42.28, which is compared with SBDD compounds. Cmp14, 22, 52, and 108 showed good binding free energy (ΔG) ­65.40, ­66.66, ­65.40, and ­65.98, respectively. Cmp17, 36, 118 and 119 also have favorable binding energies of ­45.69, ­58.63, ­54.93, and ­59.77, respectively. All the above compounds have promising free energy calculated in the HDAC8 docked complexes (Table 4). Using the compounds designed by SBDD, the physiochemical properties were calculated in Qikprop simulation. All the compounds obey the Lipinski's rules: molecular weight below 500 Da, hydrogen bond donor (less than five except for Cmp36) and acceptor (less than ten). QPlogPo/w (octanol/water partition coefficient) for all the compounds is less than five. Cmp36 alone is violating (Table 4). Total solvent accessible surface area (SASA), hydrophobic component of the SASA (FOSA) and hydrophilic component of the SASA (FISA) were analyzed for the compounds that were abiding the ranges in Qikprop physiochemical properties. Qualitative human oral absorption was predicted. Cmp14, 17, 22, and 52 have high oral absorption (score: 3). Cmp108, 118, and 119 have medium oral absorption (score: 2). However, Cmp36 alone has low oral absorption (score: 1). Polar nitrogen and oxygen van der Waals surface area (PSA) of SBDD compounds were fulfilling the limit in physicochemical calculation (Table 4). All the compounds satisfy the values of partition coefficient of octanol/gas (QPlogPoct), water/gas (QPlogPw) and brain/blood (QPlogBB). Aqueous solubility (QPlogS) and skin permeability (QPlogKp) predicted for ligands lie in the allowed the solubility and permeability range (Table 4).

Conclusion:

In this study, 120 hydroxamic acid derivative inhibitors designed using SBDD approach were screened and eight compounds were taken against HDAC8 for the docking studies. These results were compared with SAHA using docking studies, hydrogen bond and hydrophobic interactions, binding free energy and physicochemical properties of ADME/Tox. Cmp22 and 119 have good docking score in IFD. Glide energy and emodel calculations showed that Cmp17 and Cmp22 have good values. Remaining SBDD compounds showed favorable docking score, glide energy and emodel. SBDD docked complexes showed good stability from energy calculations. Hydrogen bond interactions were analyzed for SBDD compounds that have favorable interactions with active site residues. Cmp36 has nine hydrogen bond interactions at the active site residues. However, the D…A distance range is 2.78­ 2.97 Å for all the hydrogen bonds in Cmp22 with the active site residues and rest of the compounds possess some of the hydrogen bonds with distances more than 3.00 Å. G151, F152, H180, F208, M274 and Y306 residues were primarily hydrophobic and form a tunnel shape at the active site of HDAC8 [26]. IFD reveals that G151, F152, H180, F208, M274 and Y306 residues have hydrophobic interactions. Octahedral coordination is formed with Zn (II) ion in all the complexes. Binding free energy calculation has favorable scores in the SBDD complexes. However, Cmp22 has good free energy value among all the complexes. In physicochemical properties, Cmp14, 17, 22, 52, 108, 118 and 119 abide the Lipinski's rule of five and human oral absorption. All the compounds satisfied surface area calculations of SASA, FOSA, FISA, PSA and partition coefficient of QPlogPoct, QplogPw and QplogBB. All the SBDD compounds (except Cmp36) have a potent inhibitory property. Cmp22 satisfies all the in silico criteria like docking score, glide energy, glide emodel, free energy score, ADME/Tox, hydrogen bond and hydrophobic interactions. All SBDD compounds are favorable against HDAC8 and Cmp22 is found to be more potent among them, suggesting that inhibition of HDAC8 by Cmp22 could be revealed through in vivo and in vitro studies so that it could be used against various cancer therapies.