Chemical and Pharmaceutical Bulletin
Online ISSN : 1347-5223
Print ISSN : 0009-2363
ISSN-L : 0009-2363
Regular Articles
An Improved Protocol for the Virtual Screening Discovery of Novel Histone Deacetylase Inhibitors
Qiuhang SongTingting LiuYucui LiuShuyue WangCong FanLihua ZhengYongli BaoLuguo SunChunlei YuYing SunZhenbo SongGuannan WangYanxin Huang Yuxin Li
Author information
JOURNAL FREE ACCESS FULL-TEXT HTML
Supplementary material

2019 Volume 67 Issue 10 Pages 1076-1081

Details
Abstract

Histone deacetylases (HDACs) are enzymes that play a key role in structural modification and gene expression. The overexpression of HDAC is associated with cancer, and thus inhibiting the enzyme could be an efficient cancer therapy. To discover new HDAC inhibitors (HDACis), we proposed an improved protocol combining a hierarchical pharmacophore search, molecular docking, and molecular dynamic simulations. The test results showed that the improved screening protocol effectively reduced the false-positive rates of drug-like chemicals. Based on the protocol, we obtained 16 hit compounds as potential HDACis from the Life Chemicals database. Enzyme inhibition experiments showed that two of the hit chemical compounds had HDAC-inhibitory effects. In vitro assays showed that Z165155756 could selectively inhibit the proliferation of cancer cells and specifically promoted apoptosis and induced G1/S phase arrest in A2780 cells. It may have potential therapeutic effects in ovarian cancer and is worthy of further investigation.

Introduction

Histone deacetylases (HDACs) can remodel chromatin and play an important role in regulating gene expression.1) To date, based on HDAC sequence similarity to yeast orthologues Rpd3, Hdal, and Sir2, the HDAC family is divided into four classes: class I (HDAC1–3, and –8); class II (HDAC4–7, 9, and 10); class III (sirtuins 1–7); and class IV (HDAC11).2) Class I, II, and IV HDACs are classical Zn2+-dependent enzymes harboring a catalytic pocket with Zn2+-chelating compounds.3) Overexpression of HDACs in various cancers may inhibit the expression of growth-suppressive genes (p21 and p27) that play an important role in promoting cancer cell proliferation.4) HDAC inhibition has therefore become an important target for cancer therapy.5) In the expanding field of anticancer drugs, HDAC inhibitors (HDACis) are playing an increasingly important role.6) At present, vorinostat (SAHA), romidepsin (FK228), belinostat (PXD-101), panobinostat (LBH589), and chidamide are clinically used in treating cutaneous T-cell lymphomas, multiple myeloma, and peripheral T-cell lymphoma.7) However, dose-limiting toxicities including thrombocytopenia, arrhythmia, gastric side effects, nausea, and fatigue are observed with the administration of these HDACis.8) Novel HDACis with fewer side effects are urgently needed.

High-throughput screening of large chemical compound databases has been a primary source for the identification of novel lead chemical compounds.9) Virtual screening is a powerful method used for the rapid discovery of novel, original chemical candidate compounds from chemical databases like that of Life Chemicals (www.lifechemicals.com).10) The method has also been used to identify many HDACis with novel structures.11,12) Currently, pharmacophore-based studies, molecular docking, and molecular dynamic (MD) simulation approaches are used independently or in combination to evaluate effective inhibitors.12) Based on our research results, we present an improved protocol for the discovery of novel HDACis involving screening of different structure-based pharmacophore models and docking into the active site of HDAC8 to find drug-like compounds in the Life Chemicals database. Complexes of compounds with high docking scores and HDAC8 are then subjected to MD simulation, followed by analysis of the interactions of these compounds with HDAC and their free energies (△Gbind) to identify potential HDACis.13)

After the suite of screening procedures, 16 compounds with high △Gbind scores were identified. Subsequent enzyme activity experiments using an HDACi screening kit indicated that 2 of the 16 candidates had HDAC-inhibitory effects. In vitro assays indicated that Z165155756 induces apoptosis and inhibits the proliferation of several types of cancer cells, especially of A2780 cells. The improved protocol for the virtual screening discovery of novel HDACis led to the identification of the chemical candidate Z165155756, which is a promising HDACi, has potential therapeutic effects in ovarian cancer, and is worthy of further investigation.

Experimental

Pharmacophore Modeling

We used the GALAHAD program in SYBYL-X 2.0 for ligand-based pharmacophore modeling, after which 7 HDACis with different structures were selected as representative chemical compounds. Most of the parameters used the default values, with the activity column of pIC50, the exception of 100 generations, and keeping the best model of 20. At least 5 of 7 features in the pharmacophore model had to match. Finally, UNITY in SYBYL-X2.0 was used for flex searches to screen the Life Chemicals candidate database.

Molecular Docking (MD)

GOLD 5.2 software and HDAC8 (PDB ID: 1T69) were selected as the molecular docking screening method and the docking target, respectively. Except for the value of the genetic algorithm (GA) using 30, other parameters were kept as defaults. The top three scoring structures of every ligand were retained at the end of the calculation. When GOLD showed a docking conformation through a specific algorithm, a scoring function was needed to evaluate the ligand receptor affinity. The scoring function took into account the ligand and the receptor of every force, such as salt bridges, hydrogen bonds, hydrophobic stacking forces, etc. In our experiments, we used ChemPLP and ChemScore as the fitness functions for screening.

MD Simulation and Molecular Mechanics—Poisson Boltzmann Surface Area Calculations

We used the co-crystallized or the docked binding orientation structures for MD simulation, and the MD simulations of the compounds were carried out using the GROMACS software. We chose the general AMBER force field (AMBER03) and the standard AMBER force field separately for the ligands and protein (1T69). Then the protein and ligands were immersed in a truncated hexahedron box of water (TIP3P) with a margin of 10 Å along each dimension, and counterions were used to neutralize the system. Finally, all the simulations were analyzed under periodic boundary conditions, with all the parameter settings in the MD simulations similar to those reported by Wang et al.14) The binding free energies were derived using the molecular mechanics—Poisson Boltzmann surface area (MM-PBSA) method integrated in g_mmpbsa.15)

HDAC Inhibitory Effect Assay

An HDACis screening kit (APExBIO, Houston, TX, U.S.A.) was used to detect the inhibitory effects of the screened candidate chemical compounds. For positive controls or candidate inhibitors, after dissolving Trichostatin A (TSA) or candidate compounds into the appropriate solvent, diluting to 2-fold the desired test concentration with double-distilled water, and adding 50 µL of the diluted candidate inhibitor into wells, the concentration of a U.S. Food and Drug Administration-approved HDACi (trichostatin A) or candidate compounds were 20 and 100 µM, respectively. Then the chemical compounds, HDAC fluorometric substrate (with an acetylated lysine side chain), assay buffer, and HeLa nuclear extracts were placed in 96-well plates and incubated at 37°C for 25 min. Lysine developer was then added, and the mixture was incubated for another 30 min to stop the reaction. The optical density (OD) value was finally calculated as: (OD value of test HDACi/OD value of water) × 100%.

Cell Lines and Cell Culture

The human hepatocellular carcinoma cell line HepG2, human ovarian carcinoma cell line A2780, normal human liver cell line L02, and human ovarian surface epithelial cell line (HOSEpiC) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, P. R. China). All cells were cultured with Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (Biyuntian, Hangzhou, P. R. China), 100 U/mL penicillin, and 100 mg/mL streptomycin (Roche, Upper Bavaria, Germany) at 37°C under a 5% CO2 atmosphere.

3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium Bromide (MTT) Assay

The MTT assay was used to analyze cell viability after treatment with the chemical candidates. First, cells were plated into 96-well plates (1 × 104 cells/well), cultured for 24 h, and then dimethyl sulfoxide (DMSO) or 0, 1, 10, 50, or 100 µM of candidate compounds in the presence of 3% serum was added to the cultures. After 44 h of culture, 20 µL of MTT (5 mg/mL) was added to each well for an additional 4-h incubation. The formazan precipitate was dissolved with 100 µL DMSO, and finally a micro-enzyme-linked immunosorbent assay (ELISA) reader was used to measure inhibition of cell viability at 570 nm. The inhibition ratio was expressed as: (1-OD value of compounds/OD value of control) × 100%.16)

Cell Cycle Arrest and Cell Apoptosis Analysis

To quantify cell cycle arrest and cell apoptosis, A2780 cells were seeded onto 6-well plates and cultured for 12 h and then incubated with 40 µM of DMSO for 48 h or with candidate inhibitors for 24, 36, and 48 h. Then, the cells were centrifuged and stained with 5 µL fluorescein isothiocyanate (FITC)-Annexin V and 5 µL PI (Biyuntian) for 20 min or fixed in 75% alcohol at 4°C for 4 h, followed by the addition of RNase and staining with PI (Biyuntian) for 20 min. Cell cycle arrest and apoptosis rates were then measured using flow cytometry (BD FACSCant, Becton Dickinson, NK, U.S.A.) and analyzed using Modfit software (Verity Software House, Topsham, ME, U.S.A.).

Statistical Analysis

The experimental results are expressed as mean ± standard deviation (S.D.) HDAC activity was normalized to the level of DMSO control activity. All experiments were repeated at least three times, and statistical analyses of the data were performed using Student’s t-test.

Results

Pharmacophore-Based Virtual Screening

We used pharmacophore modeling methods to extract the steric and electronic features of ligand–receptor interactions and then used those features for rapid screening of the Life Chemicals database.10) Seven known HDACis from the literature were used to generate 20 pharmacophore models with GALAHAD, and the top 5 models were chosen for database screening17) (Supplementary Table 1). All models derived from more than six ligands (N_NITS ≥ 6 and FEATS ≥ 6), representing a different trade-off among the conflicting demands of maximizing steric consensus (STERICS), maximizing pharmacophore consensus (HBOND), and minimizing energy (ENERGY) (Supplementary Table 2). We used the top five models, MODEL_011, MODEL_002, MODEL_018, MODEL_014, and MODEL_007, to validate their screening abilities. The five pharmacophore models were used to screen two databases: the HDAC8 inhibitor dataset, which contains 39 active inhibitors and 1,521 unbiased decoys; and the Maximal Unbiased Benchmarking Data Sets for HDACs (MUBD-HDACs),16) which contain all the classical HDACs, including 631 HDACis and 24,609 unbiased decoys. Parameters such as the hit list (Ht), active hits (Ha), and enrichment factor (EF) were used to calculated the pharmacophore model. The results indicated that the top three models with the best EF values were MODEL_002, MODEL_011, and MODEL_018 (Supplementary Table 3). MODEL_002 and MODEL_011 consist of 6 pharmacophore features, three hydrophobes (HY4, HY5, and HY6), two hydrogen bond (HB) acceptors (AA3 and AA2), and an HB donor (DA1). MODEL_018 consists of three hydrophobes (HY5, HY6, and HY7), three HB acceptors (AA3, AA4, and AA8), and one HB donor (DA1) (Fig. 1A). All the models contain the same pharmacophoric features, but MODEL_018 has an additional HB acceptor feature. Although the chemical features of the models are similar, they differ in their 3D structures, and the differences are considered to be the major advantage of the dynamic structure-based pharmacophore models.

Fig. 1. Generated Pharmacophore Model and △Gbind Values of Candidate Chemical Compounds

(A) Pharmacophore models of MODEL_002 and MODEL_011 (length unit: angstrom). Both have three hydrophobes (HY4, HY5, and HY6), two hydrogen bond (HB) acceptors (AA3 and AA2), and an HB donor (DA1). The pharmacophore model of MODEL_018 has three hydrophobes (HY5, HY6, and HY7), three hydrogen bond (HB) acceptors (AA3, AA4, and AA8) and an HB donor (DA1). Gray spheres, hydrophobes; light gray spheres, HB acceptors; dark gray spheres, HB donors. (B, C) g_mmpbsa was used for calculating the △Gbind values of compounds and HDAC8. Light gray bars, active compounds; white gray bars, inactive compounds. * p < 0.05, Student’s t-test.

The three pharmacophore models were used to screen the Life Chemicals database consisting of 617822 small molecules. First, we used MODEL_018 with the smallest EF value of the three pharmacophore models to screen the database, and 456541 candidates were hit. Then we used MODEL_011 for screening, and 201270 candidates were hit; finally, we used MODEL_002 for screening, and 98122 candidates were hit. However, when only MODEL_002 was used to screen the database, the number of molecules was 4-fold greater than that using the hierarchical screening method. These results indicated that a high-efficiency hierarchical screening method using pharmacophore models generated by GALAHAD can be used to screen for novel HDACis.

Molecular Docking-Based Virtual Screening of Novel HDAC8 Inhibitors

In order to estimate the rationality of the 98122 candidates hit by pharmacophore-based screening and confirm whether they could bind to the active pocket of HDACs, all the candidates were docked into the active site of HDAC8 using a semi-flexible docking method reported in a previously validated protocol.17) When the GA value is 30, the average RMSD value is 1.08 Å and the docking scores derived from Chemscore and ChemPLP were used for docking. According to the data, 463 small candidates were screened out when we set the cutoff values of the ChemPLP score at >68.5 and the Chemscore at >21.4.17) Finally, we used SYBYL-X 2.0 for chemical diversity analysis, and 31 hit candidate compounds with different structures were screened out and used for further validation (Supplementary Table 4).

MD-Based Virtual Screening of Hit Compounds with HDAC8

Recently, MD simulation has been used for analyzing biomolecular interactions with the Poisson–Boltzmann surface area (MM-PBSA) method to predict the interactions of drug-like compounds with HDACs and analysis the △Gbind, and higher △Gbind meant stronger association,1619) it is an essential validation before drawing conclusions from the docking results. We collected 10 different types of HDACs from the PBD database (1T69, 1T64, 4QA3, 3F07, 3ZNR, 3C10, 4CBY, 4LXZ, 3MAX, and 4BKX) for preliminary testing in which the active sites were combined with different crystalline ligands. Moreover, 10 decoys (ZINC06430861, ZINC03413613, ZINC55476972, ZINC09112458, ZINC58202414, ZINC58381831, ZINC02933910, ZINC43736119, ZINC08254751, and ZINC10760033) were randomly selected from the MUBD database as decoys.20) The results showed that there was a significant difference between active sites and decoys (Fig. 1B), indicating that potential inhibitors of HDAC may have higher free energy. Therefore, MD simulation could be used to identify new HDACis.

MD simulation was next used to identify effective compounds after docking screening. The △Gbind values of active sites and decoys are listed in Supplementary Table 5. Compounds with △Gbind values >9.25 (average of decoys) were selected as effective candidate compounds, and 14 small molecules (Z165155756, Z46582199, Z1021075944, Z234820564, Z86187074, Z152354082, Z169686860, Z226347122, Z56765091, Z18866337, Z126203970, Z226781892, Z17378411, and Z126202522) were identified (Fig. 1C). The hit compounds were purchased from Life Chemicals for further in vitro validation.

Inhibitory Enzyme Activity Evaluation

After virtual screening, in vitro experiments were conducted to analyze the enzyme inhibitory activity of the final 14 hit chemical compounds. The results indicated that two compounds, Z165155756 and Z46582199, can significantly inhibit the enzyme activity of HDACs. After treatment with the hit candidates, the relative enzyme activities were 54% and 55%, respectively, but the other 12 small molecules did not significantly inhibit the activity against HDACs (Fig. 2A).

Fig. 2. Inhibitory Effects of the 14 Hit Chemical Candidates against HDACs and Molecular Docking Results of Two Hit Compounds

(A) Inhibitory effects of the hit HDACi. The positive control (DMSO) is shown in gray, the inhibitor control in white gray, and candidate compounds in light gray. Results are expressed as mean ± S.D. (n ≥ 3). * p < 0.05. (B, C) Docking analysis of Z165155756 and Z46582199 to HDAC. Lines, amino acid residues; spheres, the metal ion (Zn2+); white gray dotted line, hydrogen bonds.

The structures of the two lead chemical candidates did not belong to either of the four main classes of HDACis, i.e., hydroxamic acids, aliphatic acids, benzamides, and cyclic peptides.21) To investigate the binding modes of the two active candidates further, we analyzed the docking positions of the two active molecules relative to the active site of HDAC8, which showed important interactions between the chemical candidates and active site residues (Fig. 2B). Z165155756 has three functional groups: phenylmethoxy; phenyl; and quinazolinyl. The quinazolinyl group can form covalent chelate complexes with ZN2+ in the active site of HDAC enzymes and lead to reversible inhibition of their biological functions. It can also form hydrogen bond interactions with HDAC residues GLY-151 and TYR-306. Z46582199 has four main functional groups: benzimidazolyl; nitrophenyl; methylimidazole; and morpholino. The morpholino group has hydrophobic contact with HDAC residue PHE208. The two lead chemical candidates were selected for further analysis of their antitumor effects.

Antiproliferative Activity and Apoptosis-Inducing Mechanism

HDACis can selectively induce cell cycle arrest and apoptosis in many cancer cells.20) To analyze the cytotoxicity of the candidates Z165155756 and Z46582199, MTT and flow cytometry assays were designed and performed. First, HepG2 and A2780 cells were treated with four different concentrations (1, 10, 50, and 100 µM) of the two candidate compounds for 48 h, after which the IC50 values indicated that Z165155756 showed stronger inhibitory effects than Z46582199 (Table 1). The IC50 values of Z165155756 against A2780 and HepG2 cells were 34.92 ± 6.22 µM and 62.29 ± 7.54 µM, respectively. To determine whether Z165155756 is nontoxic to normal cells, we also treated HOSEpiC and L02 cells with different concentrations of the candidate compounds for 48 h, after which the IC50 values were 189.86 ± 22.07 µM and 224.39 ± 28.12 µM, respectively. The IC50 values of Z46582199 against the two cancer cell lines were greater than 300 µM, and therefore it may not exert anticancer activity in these cells. These results indicated that Z165155756 shows promising selective inhibitory activity against cancer cell viability, particularly in A2780 cells.

Table 1. Comparison of the IC50 Values of Z165155756 and Z46582199 against the A2780, HOSEpiC, HepG2 and L02 Cell Lines
ChemicalsIC50 (µM)
A2780H0SEPICHepG2L02
Z16515575634.92 ± 6.22189.866 ± 22.0762.29 ± 7.54224.39 ± 28.12
Z46582199>300

Because Z165155756 significantly inhibited the viability of A2780 cells, we also analyzed the changes after Z165155756 co-culture using flow cytometry. To determine whether the inhibition of cancer cell viability is associated with the induction of cell cycle arrest, A2780 cells were treated with Z165155756 for 0, 24, 36, and 48 h and then stained with PI. The percentage of cells in G1/S phase increased significantly by 17.95 ± 3.1%, 17.29 ± 2.6%, and 15.1 ± 4.2% at 24, 36, and 48 h, respectively (Fig. 3A). In addition to cell cycle arrest, induction of apoptosis is another key mechanism of HDACis.22) The percentage of apoptotic A2780 cells after co-culture with Z165155756 for 0, 24, 36, and 48 h was also examined and found to be 6.8 ± 1.2, 12.9 ± 2.7, 27.1 ± 5.6, and 63.9 ± 4.2%, respectively (Fig. 3B), indicating that Z165155756 induces apoptosis in a time-dependent manner. The results showed that Z165155756 markedly inhibits the viability of A2780 cells via selective induction of G1/S phase arrest and promotes apoptosis.

Fig. 3. Z165155756 Induces G1/S Arrest and Apoptosis in A2780 Cells

(A) A2780 cells were incubated with DMSO (40 µM) for 48 h or Z165155756 (40 µM) for 24, 36, or 48 h. Harvested cells were stained with PI, and the cell cycle distribution was analyzed with flow cytometry. (B) A2780 cells were incubated with DMSO (40 µM) for 48 h or Z165155756 (40 µM) for 24, 36, or and 48 h. Harvested cells were stained with Annexin V-FITC/PI, and the percentage of apoptotic cells was analyzed with flow cytometry. Three independent experiments were performed, and the results are expressed mean ± S.D.

Discussion

We improved a hierarchical virtual screening protocol through the screening of different structure-based pharmacophore models and docking to the active site of HDAC8 to identify drug-like compounds in the Life Chemicals database. Compounds with high docking scores were used in MD simulations to analyze the interactions of hit compounds with HDAC. Potential HDACis were identified based on the △Gbind values. The improved protocol appears to be an essential validation before drawing any conclusions from the pharmacophore and docking results. Based on the improved virtual screening workflow, 14 hit chemical compounds appeared to be potential HDACis. The 14 hit chemical candidate compounds were assessed using SciFinder scholar (https://scifinder.cas.org/) and the results indicated that they had not previously been assessed for their HDAC-inhibitory effects. The in vitro enzyme inhibition assay confirmed that Z165155756 and Z46582199 have significant HDAC-inhibitory effects. The two lead chemical candidate compounds have lipophilic linker parts and zinc-chelating groups, with the same hydrophobic features and HB acceptors as in the pharmacophore models. Moreover, the two candidate compounds have novel structures that do not belong to the four classes of HDACis, offering opportunities for future HDACi design.

The cytotoxicity assay indicated that, in agreement with its inhibitory effects on HDAC, Z165155756 shows promising selective inhibitory activity against cancer cells, particularly A2780 cells. Moreover, Z165155756 significantly inhibited the viability of A2780 cells via selective induction of G1/S phase arrest and promoted apoptosis. All these results indicate that the method reported here improves the HDACi screening protocol. The method was used to identify the novel structure of the promising HDACi Z165155756. Based on its effects in human ovarian carcinoma cells, Z165155756 is worthy of further investigation.

Acknowledgment

This work was supported by the Natural Science Foundation of Jilin Province, P. R. China (No. 20180101242JC), Natural Nature Science Foundation of China (81700709), Research Foundation of Jilin Provincial Science&Technology Development (20180520105JH), and Systems Biology Research on Genome and Transcriptome of Stem Cells (2017030) of Jilin Province Sunbird Regenerative Medical Engineering Co., Ltd.

Conflict of Interest

The authors declare no conflict of interest.

Supplementary Materials

The online version of this article contains supplementary materials.

References
 
© 2019 The Pharmaceutical Society of Japan
feedback
Top