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.[15] 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,16,17] 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.[18−22] 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 Zn2+-bound water molecule.[23] Another inhibitor, ACY-1083, coordinates to Zn2+ and interacts with Y745. The C=O group makes a hydrogen bond with the water molecule linked with Zn2+; 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 Zn2+, 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.[24,25] 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.[26] Molecular modeling techniques are an effective approach to understand protein–ligand interaction dynamics.[27,28] 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 Schrödinger suite were used to modify the deficiencies of the proteins[29] 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.[30,31] The ligand file preparation was carried out using Ligprep.[32] Docking was completed with the Glide tool of the Schrödinger 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.[33,34] 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.[35] 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.[36] 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/mol2 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 Nosé–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 root-mean-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.[37,38] For the binding energy investigation, Prime-MM/GBSA (molecular mechanics/generalized Born surface area) was employed.[39] This model uses a Gaussian surface rather than a van der Waals surface to represent the solvent accessible surface area.[40,41] The equations used for binding energy (ΔGbind) calculations are the following:where Ecomplex, Eprotein, and Eligand are the minimized energies for the protein–ligand complex, protein, and ligand, respectively;where Gsolv-complex, Gsolv-protein and Gsolv-ligand are the solvation energies for the protein–ligand complex, protein, and ligand, respectively; andwhere GSA-complex, GSA-protein, and GSA-ligand are the surface area energies for the protein–ligand complex, protein, and ligand, respectively.

In all docking postures, the Prime program of the Schrödinger 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.[42−44] 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 . 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 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/mol, respectively (Table ). The docking of 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 ).

Table 1

Potential HDAC6 Inhibitors of ZINC Database Docking and MM/GBSA Study

Figure 1

Molecular docking of (a) HDAC6 with Panobinostat [6918837] and (b) HDAC6 with ZINC000002845205.

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

Figure 2

RMSD profile of HDAC6 during 150 ns MD simulation (blue, HDAC6; orange, HDAC6 + Panobinostat; gray, HDAC6 + ZINC compound).

Figure 3

RMSD profile of Panobinostat and ZINC compound during 150 ns MD simulation (blue, Panobinostat; orange, ZINC compound).

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

Figure 4

RMSF study to show local changes around the protein chain (blue, HDAC6; orange, HDAC6 + Panobinostat; gray, HDAC6 + ZINC compound).

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 .

Figure 5

Compactness of the protein in terms of radius of gyration (Rg) during 150 ns MD simulation (blue, HDAC6; orange, HDAC6 + Panobinostat; gray, HDAC6 + ZINC compound).

During simulation, neither complex showed any substantial structural changes. The resulting MD trajectory was used to investigate the hydrogen-bond interactions between HDAC6 and its ligands. Figures –13 show the interactions of the H-bond formations and amino acid residues between HDAC6 and the ZINC compound throughout the 150 ns study. The top sections of Figures and 13 designate the interactions of HDAC6 with the inhibitors as a whole. The bottom sections display the interactions between specific amino acids and the ligands (Panobinostat and ZINC000002845205) in each trajectory. Dark shading indicates more than one specific contact. H-bond formations were calculated for HDAC6 in relation to the ligands under study within the 150 ns trajectory shown in Figures –13. Complex 1 showed two hydrogen bonds, wherein the first H-bond was generated by the −NH of Panobinostat with the Gly582 of the receptor with 91% residency. The second H-bond was formed by the −N of His573 with 50% residency. The amino acids Phe583 and Asn643 were involved in hydrophobic interactions. Hydrogen-bond interactions were also identified for Complex 2 after MD simulation. The first H-bond was generated by the −NH of ZINC000002845205 with the Gly582 of the receptor with 96% residency. The −N of His573 formed the second hydrogen bond with 96% residency. Hydrophobic interactions were present between amino acids Phe583 and Asn643 with ZINC000002845205. An H-bond was generated by Panobinostat and ZINC000002845205 with Gly582, and His573 substantiated the docking result. The average number of H-bond formations between HDAC6 with Panobinostat and with ZINC000002845205 were 2.12 and 2.07 respectively. The analysis of the HDAC6 Panobinostat MD trajectory reveals that the interactions detected in the structure are preserved all over the MD simulation. In conclusion, the MD simulations suggest that the Panobinostat makes hydrogen bonds with the Gly582 and His573 in HDAC6.

Figure 6

H-bond formation between HDAC6 + Panobinostat during 150 ns MD simulation study.

Figure 13

Top block shows specific contacts, and the bottom block displays HDAC6 + ZINC compound interactions during 150 ns simulation (dark orange shade represents a stable interaction).

Figure 11

Top block shows specific contacts, and the bottom block displays HDAC6 + Panobinostat interactions during 150 ns simulation (dark orange shade represents a stable interaction).

Figure 10

HDAC6 + Panobinostat interaction during 150 ns simulation (water bridges, blue; H-bond, green; hydrophobic contacts, violet).

A 2D schematic of ligand protein contacts (HDAC6 + Panobinostat) during 150 ns MD simulation study.

HDAC6 + ZINC compound interactions during 150 ns simulation (water bridges, blue; H-bond, green; hydrophobic contacts, violet).

IC50 values of 3–61 nM were reported for Panobinostat with HDAC1–HDAC9, excepting for a slightly higher IC50 value of 248 nM for HDAC8.[45,46] 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 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 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.

Figure 9

A 2D schematic of protein ligand contacts (HDAC6 + ZINC compound) during 150 ns MD simulation study.

Figure 8

H-bond formation between HDAC6 + ZINC compound during 150 ns MD simulation study.

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. Docking and MD simulations support the computation of the free energy of the HDAC6 + ZINC000002845205 binding by using theoretical calculations.[47] In free energy calculations, docking and scoring via MM/PBSA [molecular mechanics (MM) with Poisson–Boltzmann (PB) and surface area solvation] lacks accuracy but is fast and can classify ligands into binders and nonbinders.[48] MM/GBSA, which is an alternative method, is more statistically accurate than the docking and scoring functions.[49−52]where Gbind is the binding free energy, H is enthalpy, T is temperature, and S is entropy.where EMM is the gas phase energy and Gsol is the solvation free energy.where Einternal is the bond, angle, and dihedral energies; Eelectrostatic is the electrostatic energy; and Evdw is the van der Waals energy.where Gsol is the solvation energy, GPB/GB is the nonelectrostatic solvation energy for polar contributions, and GSA is the electrostatic solvation energy.where Gnonpolar 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 ΔGPB/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 . It was evident that ZINC000002845205 showed promising results compared with known inhibitor Panobinostat (Table ).

Table 2

Binding Energy (kcal/mol) of HDAC6 with Panobinostat and ZINC000002845205

	protein + Panobinostat	protein + ZINC000002845205
electrostatic	–41.5497	–46.612
H-bond	–1.64344	–1.34565
van der Waals energy	–39.0996	–30.8053
lipophilic energy	–17.3848	–18.7312
π–π packing correction	–6.53995	–3.23679
solv GB	46.63933	49.99878
ΔGbind	–54.115	–47.2366

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 PCaco-2 = 1639.51 nm/sec[53] and a PMDCK = 844.17 nm/sec.[54] The predicted in vitro skin permeability of ZINC000002845205 was −1.94 cm/hour, which is considered to be within an acceptable range.[55] 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.[56]

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 ΔGbind reported from MM/GBSA calculations indicated that Complex 2 has a comparable binding energy to Complex 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.