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][5][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][10][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 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.

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 OPLS-AA 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).

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).
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.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 Å.  where, vdW-van der Waal energy; Coul -Coulomb energy; Lipo -lipophilic contact term; HBond -hydrogen-bonding term; Metal -metal-binding term; BuryP -penalty for buried polar groups; RotB -penalty for freezing rotatable bonds; Site -polar interactions at the active site; and the coefficients of vdW and Coul are: a =0.065, b = 0.130.Glide energy is combination of the posed ligand and receptor.Glide Emodel is a specific combination of Docking Score , CvdW (CvdW = Coul + vdW is the non-bonded interaction energy between the ligand and the receptor) and the internal torsional energy of the ligand conformer.∆G Binding free energy calculation.†Estimated number of hydrogen bonds that would be accepted by the solute from water molecules in an aqueous solution.Values are averages taken over a number of configurations, so they can be non-integer.$ Number of Lipinski's rule of five violations.SASA: Total solvent accessible surface area, FOSA: Hydrophobic component of the SASA (saturated carbon and attached hydrogen), FISA: Hydrophilic component of the SASA (SASA on N, O, and H on heteroatoms), PSA: Van der Waals surface area of polar nitrogen and oxygen atoms.QPlogPoct, QPlogPw, QPlogS and QPlogBB were predicted partition coefficient of octanol/gas, water/gas, aqueous solubility and brain/blood, respectively.QPlogKp was predicted skin permeability (log Kp).
[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].

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

Table 1 :
[26], 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, 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.Induced fit docking score, glide energy and glide emodel calculations Compounds Docking score Glide energy (kcal/mol) Glide emodel (kcal/mol)

Table 3 :
Coordination distance between coordinating residues/atoms with Zn ion

Table 4 :
Interactions, binding free energy and physicochemical properties of docked complexes