Discovery of Histone Deacetylase Inhibitor Using Molecular Modeling and Free Energy Calculations

The histone acetylation–deacetylation at lysine regulates the functions of many cellular proteins. An increased expression of HDAC6 can cause an increased amount of deacetylated histones, which leads to an inhibition of gene expression and has been associated with cancer cell proliferation. The present study screened the ZINC database to find novel HDAC6 inhibitors using virtual high-throughput screening techniques. The docking score, free energy, and binding pattern of the complexes were used to select a best ligand for further study. Molecular dynamic simulations, binding interactions, and the stability of docked conformations were investigated. Several parameters that determine protein–ligand interactions, such as root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), and binding pattern, were observed. Hydrogen bonds were observed at His 573 and Gly 582 after a 150 ns simulation with identified compound ZINC000002845205, and they were similar to known inhibitor Panobinostat. The molecular mechanics with generalised Born and surface area solvation (MM/GBSA) free energy was comparable to known inhibitor Panobinostat. ZINC000002845205 qualifies drug-likeness according to Lipinski’s rule-of-five, rule-of-three, and the World Drug Index (WDI)-like rule, but there is one violation in the lead-like rule.


■ INTRODUCTION
Lysine acetylation−deacetylation processes regulate the activity of biological proteins. The important classes of proteins involved in this cycle are lysine acetyltransferases, a recognizing protein, and histone deacetylases. Post-translational changes include reversible acetylation and deacetylation of a histone by histone acetylase and deacetylase enzymes, which add or remove acetyl groups from specific lysine amino acids. 1−6 Alterations in the interaction between histones and DNA regulate the compactness of chromatin and gene expression. 7,8 The DNA histone deacetylation causes condensation of chromatin that results in gene expression silencing and repression of gene transcription factors. 9,10 The human histone deacetylase 6 (HDAC6) has two catalytic parts comprising a ubiquitin-binding domain and a dynein-binding domain. 11−13 The active site of the HDAC6 protein is mainly formed by H463, P464, F583, and L712. 14 An increased expression of HDAC6 can cause an increase in deacetylated histones, which leads to an inhibition of gene expression and has been associated with cancer cell proliferation. 12,15,16 HDAC6 is the cytosolic tubulin deacetylase protein, and inhibition of HDAC6 can inhibit microtubule functions to cause cell cycle apoptosis. Dysregulation of HDAC6 functions, that is, abnormal expressions of HDAC6, has been associated with different kind of malignancies. 17−21 All these factors make HDAC6 an important target for anticancer therapy and other disorders. HDAC6 inhibition can cause increased acetylation on the histone lysine residues and activation of target genes that are selectively repressed in tumors, which primes cell cycle arrest. Romidepsin, Vorinostat, Belinostat, and Panobinostat are permitted HDAC6 inhibitors. HDAC6 inhibitors are in use for different therapeutic indications like cutaneous T-cell lymphoma (Vorinostat), peripheral T-cell lymphoma (Romidepsin), refractory peripheral T-cell lymphoma (Belinostat), and multiple myeloma (Panobinostat). The crystal structure of the zCD2, a surrogate for the human enzyme, showed that the N-hydroxy-4-[(N(2-hydroxyethyl)-2-phenylacetamido) methyl-benzamide] (HPB), a known HDAC6 inhibitor, binds to D612, H614, D705, and a water molecule at the catalytic site. H573 and H574 establishes hydrogen bonds with the Zn 2+bound water molecule. 22 Another inhibitor, ACY-1083, coordinates to Zn 2+ and interacts with Y745. The CO group makes a hydrogen bond with the water molecule linked with Zn 2+ ; amino acids H573 and H574 also make hydrogen bonds with this water molecule. The aromatic ring of the aminopyrimidine linker is compressed between F583 and F643, which is the same as observed with inhibitor HPB. The phenyl group of the inhibitor creates van der Waals interactions with amino acids P464 and F583. The hydroxamate part of Ricolinostat coordinates to Zn 2+ , and the side chain of Y745 donates a hydrogen bond to the hydroxamate CO group; H573 donates a hydrogen bond to the hydroxamate N−O − group, and H574 accepts a hydrogen bond from the hydroxamate NH group. 23,24 Pan-HDAC6 inhibitors, trichostatin A (TSA), H573, and H574, establish hydrogen bonds with the N−O group, while Y745 provides a hydrogen bond to the hydroxamate CO. 25 Molecular modeling techniques are an effective approach to understand protein−ligand interaction dynamics. 26,27 The present study aims to screen a novel potential HDAC6 inhibitor from ZINC databases for drug discovery and development.

■ METHODS
HDAC6 (PDB ID: 5EF8) and about 832 686 compounds of molecular weight from 300 to 350 Daand lipid solubility from 2.0 to 3.0 were retrieved from the Protein Data Bank (PDB) and ZINC ligand database, respectively. The Epik module and Protein Preparation Wizard tool of the Schrodinger suite were used to modify the deficiencies of the proteins 28 and protonation at biological pH. Protein minimization and hydrogen bond optimization were done by eliminating water molecules that were greater than 3 Å away except one water molecule in the active site. The OPLS3 (Optimized Kanhesia for Liquid Simulations) force field was used for restrained molecular minimization to make the structures relax and minimize steric clashes. 29,30 The ligand file preparation was carried out using Ligprep. 31 Docking was completed with the Glide tool of the Schrodinger suite. A cocrystallized ligand was identified, and a grid was generated around the ligand in HDAC6. Ligands under study were then docked with HDAC6. The docking was performed by exploiting default parameters for an extra precision (XP) study. 32,33 All compounds underwent molecular docking, and then the docking score, free energy, and visual inspection of the binding pattern of the complexes were used to select the best ligand for further study.
Molecular Dynamics (MD) Simulations. MD simulations of the identified protein ligand complex were done with the Desmond module of the Schrodinger suite. 34 The OPLS3 force field was used to solvate the docked HDAC6−ZINC compund unit, employing the SPC (simple point charge) aqueous solvation of the orthorhombic solvent box and a solvent buffer spreading 10 Å away from the protein. 35 Then, the system was neutralized by using appropriate counterions. Prior to production runs, the six-step relaxation protocol was applied. With and without a restriction of 50 kcal/mol 2 on the solute atoms, the first two stages of steepest descent minimization (2000 steps) were employed. The final four short MD simulations were 12 ps in an NVT ensemble at 10 K, an NPT ensemble at 10 K with the same restraint for 12 ps, an NPT ensemble at 300 K with restraint for 12 ps, and a 24 ps NPT ensemble at 300 K without any restrictions. With a relaxation time of 2 ps, the Nose−Hoover chain thermostat and the isotropic Martyna−Tobias−Klein barostat were used. Short-range interactions were considered with a cutoff of 9 Å. The smooth particle mesh Ewald method (PME) was used to examine long-range Coulombic interactions. For nonbonded interactions, the r-RESPA integrator was used, with short-range forces updated every step and long-range forces updated every three steps. For the protein−ligand complex, 150 ns MD simulations were run, and trajectories were generated every 150 ps with an energy recording interval of 1.2 ps. The rootmean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) of the structure of the protein−ligand complex were investigated using the OPLS3 (Optimized Kanhesia for Liquid Simulations) force fields with regard to a 150 ns simulation. 36,37 For the binding energy investigation, Prime-MM/GBSA (molecular mechanics/generalized Born surface area) was employed. 38 This model uses a Gaussian surface rather than a van der Waals surface to represent the solvent accessible surface area. 39,40 The equations used for binding energy (ΔG bind ) calculations are the following: where E complex , E protein , and E ligand are the minimized energies for the protein−ligand complex, protein, and ligand, respectively; where G solv-complex , G solv-protein and G solv-ligand are the solvation energies for the protein−ligand complex, protein, and ligand, respectively; and where G SA-complex , G SA-protein , and G SA-ligand are the surface area energies for the protein−ligand complex, protein, and ligand, respectively. In all docking postures, the Prime program of the Schrodinger software suite was used to calculate the MM/ GBSA. The directionality of the hydrogen-bond and π-stacking interactions was demonstrated using the variable dielectric solvent model VSGB 2.0. The complex was minimized for molecular mechanics by letting residues within 5.0 Å of the ligand to relax while keeping the rest of the structure fixed. 41−43 Drug-likeness of the most active compound was calculated using the QikProp & PreADMET programs.

■ RESULTS & DISCUSSION
The intention of this study is to recognize novel probable HDAC6 inhibitors from the ZINC database. This was achieved through a molecular modeling study of ZINC database compounds with Glide software. Redocking of a cocrystallized inhibitor with an acceptable RMSD value was used to validate the docking technique. The native and redocked poses demonstrated the same binding pattern. The docking approach is valid and reliable for predicting potential inhibitors in our compound database, as demonstrated by this run.
The screening was conducted on the basis of XP molecular docking results for 832 686 compounds at the HDAC6 active site, and the results are shown in Table 1. ZINC000002845205 displayed a good contact at the binding site by creating hydrogen bonds and hydrophobic interactions with HDAC6 and a good MM/GBSA binding energy. The MM/GBSA binding energy estimates the free energy change during ligand target binding. The hydrogen bonding occurred at Gly582 and His573. Hydrophobic bond formations occurred with amino acids Phe583 and Asn643. Figure 1 displays interactions with Panobinostat (known inhibitor) and ZINC000002845205. The manner of interaction for amino acids in both the cases were similar. The postdocking Glide score of HDAC6 with ZINC000002845205 was found to be −9.70, compared with −11.13 for HDAC6 with Panobinostat, while the MM/GBSA binding free energy was found to be −52.03 and −56.08 kcal/ m o l , r e s p e c t i v e l y ( T a b l e 1 ) . T h e d o c k i n g o f ZINC000002845205 showed two hydrogen bonds with Gly582 and His573 and hydrophobic interaction with Phe583 and Asn643 at the active site, which is similar to known inhibitor Panobinostat (Figure 1). This docking study only provides an initial understanding of the interactions and binding mechanism, but additional factors such as the solvent impact, protein flexibility, and structural stability must be considered before reaching a final judgment. ZINC000002845205 was selected for additional molecular dynamic simulation study with HDAC6 on the basis of its docking score, free energy, and polar and nonbonded interactions.
MD simulations were run for 150 ns to study the binding mode and to assess the stability of the docked conformation of the most active compound, i.e., ZINC000002845205. Molecular dynamics simulations explored the dynamic interactions, structural stability, and binding site adjustments to the docked Panobinostat and ZINC000002845205 during the 150 ns simulation. HDAC6 complexes were investigated to check the fluctuations and stability of the interactions during the 150 ns simulation using the backbone atoms of proteins and ligands (Panobinostat and ZINC compound). Figure 2 shows the RMSD for Complex I (HDAC6 + Panobinostat) and Complex 2 (HDAC6 + ZINC000002845205). For Complex 1, the average RMSD for the HDAC6 and Panobinostat are 1.13 and 1.47 Å, respectively. The RMSD of the HDAC6 and ZINC000002845205 of Complex 2 are 1.21 and 0.65 Å. Complex 1 reached equilibrium after 5 ns and remained stable afterward. The average RMSD of the ZINC000002845205 compound is comparable to the known inhibitor Panobinostat (Figure 3).
The root-mean-squared fluctuations (RMSF) of the amino acids surrounding the ligand were estimated during MD simulation to check the stability of the active site. The HDAC6 structures of the two complexes have similar RMSF and dynamic characteristic trends. The RMSF values of the residues around the ligand, Gly582 and His573, were analyzed in relation to the initial structures. The RMSF for each residue surrounding the ligands is less than 0.53 Å in all of the complexes. RMSDs of around 1.21 Å, with RMSF fluctuations of less than 0.53 Å, showed the stability of the interactions (Figure 4).
The radius of gyration interactions between HDAC6 and the ligands were also studied. This measure indicates the molecule's shape at any given point in time. It was computed when the Panobinostat and ZINC000002845205 were bound in order to examine structural alterations in HDAC6. The graph of the radius of gyration in simulation time for the Panobinostat and ZINC000002845205 complex is given in Figure 5.
During simulation, neither complex showed any substantial structural changes. The resulting MD trajectory was used to investigate the hydrogen-bond interactions between HDAC6            44,45 The HDAC6 + ZINC000002845205 complex represented a similar kind of binding pattern. Evaluation of the initial and final structures displayed a stable docked conformation during the entire simulation of 150 ns (Figures 9 and 13). The hydrogen-bond construction with amino acids Gly582 and His573 were found to be constant with the ligand. The residues Gly582 and His573 were critical for the inhibitors' binding, according to the entire analysis of the evaluated trajectories. The variations of distance between the backbone atoms of the important amino acids at the active site region with ZINC000002845205 were estimated and plotted (Figures 6 and 8). It was evident that hydrogen-bond formation is a crucial event, as indicated by the stable interaction in the docking and molecular dynamics results.
Study of the free energy gives an insight on interaction, folding pattern, and other mechanisms. The components included in free energy determination are various bonded and nonbonded interactions, such as ionic, hydrogen, electrostatic, and van der Waals, and polarization of the interacting groups.
where G bind is the binding free energy, H is enthalpy, T is temperature, and S is entropy.
where E MM is the gas phase energy and G sol is the solvation free energy.
where E internal is the bond, angle, and dihedral energies; E electrostatic is the electrostatic energy; and E vdw is the van der Waals energy.
where G sol is the solvation energy, G PB/GB is the nonelectrostatic solvation energy for polar contributions, and G SA is the electrostatic solvation energy.
where G nonpolar is the nonelectrostatic solvation energy for nonpolar contributions, SASA is the solvent-accessible surface area, γ is a correlation coefficient, and b is a fitting parameter. Polar contributions to the equation are represented by ΔG PB/GB for the GBSA or PBSA models, while SASA represents nonpolar energy contributions. The calculation of −TΔS was done from snapshots retrieved during the MD simulations. The MM/GBSA binding energy assessment was performed for each Complex. The binding energies of HDAC6 + Panobinostat and HDAC6 + ZINC000002845205 are shown in Table 2. It was evident that ZINC000002845205 showed promising results compared with known inhibitor Panobinostat ( Table 2).
The absorption, distribution, metabolism, and excretion (ADME) of drugs are important components of pharmacokinetics. According to statistics, many drug candidates miss clinical trials because of ADME issues; hence, it is important to research ADME qualities before selecting compounds as therapeutic candidates. Predicted qualitative human oral absorption was 3 (high), and percent human oral absorption on a scale of 0−100% was found to be 96% for ZINC000002845205, which falls under the category of good absorption (QikProp). Human intestinal absorption (HIA) was predicted as good at 100% for ZINC000002845205. Caco-2 and MDCK cell models are considered as reliable in vitro models to predict oral drug absorption. Caco-2 cells are derived from human colon adenocarcinoma and have multiple drug transport pathways through the intestinal epithelium. MDCK cell refers to Martin−Darby canine kidney cells. The compound was also of a high permeability group with a P Caco-2 = 1639.51 nm/sec 52 and a P MDCK = 844.17 nm/sec. 53 The predicted in vitro skin permeability of ZINC000002845205 was −1.94 cm/hour, which is considered to be within an acceptable range. 54 The predicted Blood−Brain Barrier (BBB) penetration, i.e., logBB, of ZINC000002845205 is also within an acceptable range (−0.70). ZINC000002845205 qualifies drug-likeness per Lipinski's rule-of-five, rule-of-three, and the World Drug Index (WDI)-like rule, but with one violation observed in the lead-like rule. 55

■ CONCLUSION
In the present study, the ZINC database was screened to find novel HDAC6 inhibitors using virtual high-throughput screening techniques. Many observations were noted during this process. Gly582 and His573 were recognized as important residues by the molecular modeling of HDAC6 with a known inhibitor, Panobinostat, and with ZINC database compound ZINC000002845205. The binding modes of Panobinostat and ZINC000002845205 with HDAC6 were comparable. The RMSDs of Complex 1, HDAC6 and Panobinostat, were 1.13 and 1.47 Å, respectively. The RMSDs of Complex 2, HDAC6 and ZINC000002845205, were 1.21 and 0.65 Å, respectively. The RMSFs at Gly582 and His573 were lower than 0.6 Å, which specifies the stability of the binding pocket throughout the MD simulations. The average number of H-bond formations of Panobinostat and ZINC000002845205 with HDAC6 were 2.12 and 2.07, respectively. The ΔG bind reported from MM/GBSA calculations indicated that Complex 2 has a c o m p a r a b l e b i n d i n g e n e r g y t o C o m p l e x 1 . ZINC000002845205 qualifies the drug-likeness property according to Lipinski's rule-of-five, WDI-like rule, and the rule-of-three, but with one violation observed in the lead-like rule. These observations give an insight about the potential of the ZINC000002845205 compound for further drug discovery and development process as an HDAC6 inhibitor.