Design, virtual screening, molecular docking and molecular dynamics studies of novel urushiol derivatives as potential HDAC2 selective inhibitors

Gene 637 (2017) 63–71

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Design, virtual screening, molecular docking and molecular dynamics studies of novel urushiol derivatives as potential HDAC2 selective inhibitors

MARK

Hao Zhoua,b,⁎, Chengzhang Wanga,b,⁎, Jianzhong Yea, Hongxia Chena, Ran Taoa a b

Institute of Chemical Industry of Forest Products, CAF, Nanjing, Jiangsu 210042, China Key Lab of Biomass Energy and Material, Nanjing 210042, Jiangsu, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Urushiol derivatives Molecular docking Molecular dynamics HDAC2 selective inhibitor

Three series of novel urushiol derivatives were designed by introducing a hydroxamic acid moiety into the tail of an alkyl side chain and substituents with differing electronic properties or steric bulk onto the benzene ring and alkyl side chain. The binding affinity toward HDAC2 of the compounds was screened by Glide docking. The best scoring compounds were processed further with molecular docking, MD simulations and binding free energy studies to analyze the binding modes and mechanisms. Six compounds, 21, 23, 10, 19, 9 and 30, gave Glide scores of − 7.9 to −8.5, which revealed that introducing eF, eCl, triazole, benzamido, formamido, hydroxyl or nitro substituents onto the benzene ring could increase binding affinity significantly. Molecular docking studies revealed that zinc ion coordination, hydrogen bonding and hydrophobic interactions contributed to the high calculated binding affinities of these compounds toward HDAC2 and that His145, His146, Gly154, Glu103, His183, Asp104, Tyr308 and Phe155 contributed favorably to the binding. MD simulations and binding free energy studies showed that all complexes possessed good stability as characterized by low RMSDs; low RMSFs of residues, moderate hydrogen bonding and zinc ion coordination; and low values of binding free energies. van der Waals and electrostatic interactions provided major contributions to the stability of these complexes. These results show the promising potential of urushiol derivatives as potent HDAC2 binding lead compounds.

1. Introduction Histone deacetylases (HDACs) are a group of zinc metalloenzymes that regulate chromatin remodeling and gene transcription by catalyzing the removal of an acetyl group moiety from the Ɛ-amino groups of lysine residues on the amino terminal tails of the core histones (Ning et al., 2013). HDACs play a pivotal role in the regulation of gene expression, cell growth, and proliferation (Hildmann et al., 2007). Overexpression of HDACs has been linked to the development of cancers in humans (Drummond et al., 2005). Thus, HDAC has become an important target enzyme for anticancer therapies. To date, 18 HDACs isoforms have been found in humans. Based on their homology to yeast HDACs, they can be divided into four classes: class I (HDAC1–3, and 8), class II (HDAC4–7, 9 and 10), class III (SIRT 1–7), and class IV (HDAC11) (Johnstone, 2002). Classes I, II, and IV are Zn2 +-dependent metallohydrolases, whereas class III enzymes are NAD+-dependent Sir2-like deacetylases (de Ruijter et al., 2003). Of note, HDAC2 belongs to the class I isoforms, which is highly conserved. It is overexpressed in most solid and hematological tumors and is highly

correlating with a worse prognosis, but it is not found in resting endothelial cells and normal organs (Rikimaru et al., 2007). Therefore, selective targeting of the HDAC2 isoform by directly inhibiting its functions has become a major focus in cancer chemotherapy (Rikimaru et al., 2007). HDAC inhibitors (HDACIs) have been recognized as a new class of anticancer agents. HDACIs enhance the level of histone acetylation and induce cell cycle arrest, differentiation and apoptosis (Hess-Stumpp, 2005). Many HDACIs have been developed and exhibit therapeutic efficacy in both preclinical and clinical trials. Trichostatin A (TSA) is an effective natural hydroxamate HDACI, and suberoylanilide hydroxamic acid (Vorinostat, SAHA) is the first HDACI approved by the FDA for the treatment of advanced cutaneous T-cell lymphoma (CTCL) in 2006. Additionally, many other small molecule agents, such as PXD-101, LBH589, MM-232 and CUDC-101, are in different phases of clinical trials for the treatment of different cancers (Fig. 1) (Abdel-Atty et al., 2014). Most of these HDACIs are, however, pan-inhibitors, and they show side effects that might limit their further clinical application, such as diarrhea, electrolyte changes, and cardiac arrhythmias, in clinical trials

Abbreviations: HDAC, histone deacetylases; HDACIs, histone deacetylases inhibitors; HDLP, histone deacetylases - like protein; MDS, molecular dynamics simulation; MD, molecular dynamics; RMSD, root mean square deviations; RMSF, root mean squared fluctuation ⁎ Corresponding authors at: Institute of Chemical Industry of Forest Products, CAF, Nanjing, Jiangsu 210042, China. E-mail addresses: [email protected] (H. Zhou), [email protected] (C. Wang). http://dx.doi.org/10.1016/j.gene.2017.09.034 Received 25 April 2017; Received in revised form 18 August 2017; Accepted 19 September 2017 Available online 20 September 2017 0378-1119/ © 2017 Published by Elsevier B.V.

Gene 637 (2017) 63–71

H. Zhou et al.

Fig. 1. Structures of some hydroxamate HDACIs and urushiol.

activity (Ma et al., 2012). As a result, the usefulness of urushiol as a potential anticancer agent has been limited. It is thought that the catechol group of urushiol is the main cause of allergic reactions and polymerization (Xia et al., 2004); therefore, many studies on chemical modifications of this moiety have been carried out to prepare structurally stable urushiol derivatives with no sensitization properties. In a previous study, we reported the replacement of the catechol group with a methylene acetal fragment to the synthesis of urushiol derivatives. These displayed excellent cytotoxicity activities against many types of tumor cells (Wang et al., 2015). Furthermore, some literature reports have shown that the catechol fragment can be changed to a phenylmethyl ether, methylsilyl ether or acetic ester structure to prepare urushiol derivatives, but their antitumor activities were not studied (Markiewitz and Dawson, 1965; Draper et al., 2002; Yamauchi et al., 1980). A recent investigation has shown that urushiol possesses significant HDACI activity (Ryckewaert et al., 2014). The structure of urushiol is similar to that of SAHA (Fig. 1), suggesting urushiol to be a promising lead compound for development of novel HDAC inhibitors. Compared to the classical structure of HDACI, urushiol has the linker and surface recognition (CAP) groups but has no metal binding (ZBG) group; therefore, it is necessary to introduce this key group of HDAC inhibition into its structure. Until now, there has been no report about the design and activity study of urushiol derivatives as potential HDACIs. Molecular docking and molecular dynamics have proven their pivotal role in the design and virtual screening of new bioactive molecules for drug discovery. Molecular docking has helped in designing and synthesizing inhibitors against enzymes with high efficiency and specificity; molecular dynamics simulation (MDS) provides virtual

(Yao et al., 2014). Therefore, development of isoform- and class-selective HDACIs to reduce undesirable side effects from off-target activity is highly in demand (Yao et al., 2014). X-ray crystallographic analysis of HDAC-like protein (HDLP) bound to SAHA revealed the active site of HDAC exists in a hydrophobic pocket, which comprises of a tube-like 11 Å internal channel, a Zn2 + center at the bottom of the internal channel, and a 14 Å internal cavity located near this channel (Wang, 2009). It is also known that the structural requirements of a typical HDAC inhibitor should include three distinct domains: metal binding (ZBG), linker and surface recognition (CAP) domains (Cai et al., 2015). The metal binding domain coordinates Zn2 + in the bottom of the active pocket of HDACs, and it appears to be a precondition for HDAC-inhibitory activity. Hydrophobicity is important for the linker domain. The surface recognition domain is essential for recognizing and binding to the rim of active pocket of enzymes (Huang et al., 2012b). All of the above provides useful structural information for the design and discovery of novel potent HDACIs. Urushiol is a natural compound that exists in the sap of the lacquer tree (Rhus verniciflua Stokes, Anacardiaceae), which is widely cultivated in northeastern Asian countries, including China, Korea, and Japan (Honda et al., 2008). Urushiol consists of o-dihydroxybenzene (catechol) coupled to a unsaturated alkyl side chain of 15 or 17 carbons (Honda et al., 2008). Studies show that urushiol has potent anticancer activities, showing excellent inhibitory effects on many kinds of tumor cells by inducing cell cycle arrest and apoptosis, inhibiting nuclear transcription factors and inducing angiogenesis of tumor cells (Hong et al., 1999). However, urushiol can cause hypersensitivity reactions, and it is prone to oxidation and polymerization, negatively affecting its

64

Gene 637 (2017) 63–71

H. Zhou et al.

atomic charges were calculated using MMFFs force field. The Glide module (Schrödinger, LLC, New York, NY, USA) was used to generate a receptor grid file that determine the position and size of the active site using the original ligand in the downloaded crystal structure to define the binding site and to simulate the binding mode of candidate molecular ligands to HDAC2. Glide docking was performed using the extraprecision (XP) mode with default protocols, and the ligands were treated as flexible. All the docked complexes were scored by the Glide scoring function, and the complexes with higher Glide-score were selected for further analysis of binding mode and critical interactions with the key amino acids.

realistic insights into mechanisms and efficacy of even multi-subunit enzymes inhibitors (Wang et al., 2013). Recently, the combination of docking and molecular dynamics (MD) has been successfully employed in rational HDACI design to improve initial docking results and to provide an understanding of inhibition mechanisms (Huang et al., 2012b); however, these computational methods have not been used in the structure optimization and design of urushiol derivatives for the development of new HDAC inhibitors. Considering the anticancer activities of urushiol and its structural similarities to HDACIs, we sought to develop some novel urushiol derivatives as potential HDAC2-selective inhibitors in order to search for more valuable candidates for the antitumor therapy. In the present paper, we present the design of a series of urushiol-based derivatives with a metal binding group at the tail of alkyl chain and different pharmacophores to the benzene ring or alkyl chain. The binding affinity and selectivity to the pocket of HDAC2 of the designed urushiol derivatives were evaluated by scoring functions, and those urushiol derivatives having relatively high scores were chosen for molecular docking and molecular dynamics simulations to further explore the binding modes and interactions of these compounds at the atomic level.

2.3. Molecular dynamics simulation The complexes of HDAC2 bound to ligands 21, 23, 10, 19, 9 and 30 were selected as the initial conformations for MD simulation. Each complex was then neutralized by adding sodium/chloride counter ions and solvated in a cuboid box of TIP3P water molecules with solvent layers 10 Å between the box edges and solute surface. Then, the partial atomic charges of the ligand and protein in complex were calculated by the restrained electrostatic potential (RESP) charge from the calculation with Gaussian09 package at the HF/6–31 g* level. The GAFF and AMBER14SB force fields were applied. All MD simulations were performed using AMBER 14. The SHAKE algorithm was used to restrict all covalent bonds involving hydrogen atoms with a time step of 2 fs. The particle mesh Ewald (PME) method was employed to treat long-range electrostatic interactions. For each solvated system, two steps of minimization were performed before a heating step. The first 4000 cycles of minimization were performed with all heavy atoms restrained with 50 kcal/(mol Å2), whereas solvent molecules and hydrogen atoms were unrestrained. Then, non-restrained minimization was carried out involving 2000 cycles of steepest descent minimization and 2000 cycles of conjugated gradient minimization. Afterwards, the whole system was first heated from 0 K to 300 K in 50 ps using Langevin dynamics at a constant volume, and then equilibrated for 400 ps at a constant pressure of 1 atm. A weak constraint of 10 kcal/(mol Å2) was used to restrain all the heavy atoms during the heating steps. Periodic boundary dynamics simulations were carried out for the whole system with an NPT (constant composition, pressure, and temperature) ensemble at 1 atm and 300 K in the production step. The production phase of each system was simulated for 20 ns.

2. Material and methods 2.1. Design of new urushiol derivatives Based on the structural characteristics of the FDA-approved hydroxamic acid HDACI SAHA and taking urushiol and its triene alkyl as the lead compound, we designed three series of novel urushiol derivatives (Fig. 2) based on three attributes. First, we introduced the hydroxamic acid moiety into the tail of the alkyl side chain of urushiol to obtain molecules possessing the zinc binding domain; meanwhile, to block the oxidative polymerization of urushiol, the ortho phenolic hydroxyl groups of urushiol were replaced by methylene acetal, two phenylmethyl ethers, two methyl ethers and acetate ester groups. This afforded the first series, compounds 1–4. Second, we introduced hydroxy, carbonyl, amino and methyl ether groups into the alkyl side chain of urushiol to afford the second series, compounds 5–8. Third, to investigate the effects of a variety of substituents on the benzene ring on the HDACI activity of urushiol, we introduced electron-donating or -withdrawing substituents, or groups with differing steric bulk into the benzene ring to afford the third series, compounds 9–30. The benzene ring substituents included halogens (eF, Cl, Br), fluorine-containing groups (trifluoromethyl, trifluoromethoxy), nitrogen-containing heterocyclic groups (thiazole, oxazole, pyrazole, pyrimidine, thiadiazole, triazole), sulfonamido groups (methylsulfonyl amino, phenylsulfonamido), acylamino groups (formamido, benzamido), and methyl, methoxyl, benzyloxy, amino, dimethylamino, hydroxy and nitro groups.

2.4. Binding energy calculations The MM/PBSA method in the AmberTools suite was used to calculate the binding energies between the ligands and the receptors. For each complex, one trajectory was generated through MD simulations from the initial docked compound pose. The trajectories at the equilibrium state were selected to calculate the representative binding energy. The free energy of binding, ΔGbind, was calculated using Eq. (1) from the free energy of the receptor-ligand complex (Gcpx) with respect to the unbound receptor (Grec) and ligand (Glig):

2.2. Virtual screening and molecular docking studies All 30 of the designed urushiol derivatives were screened by molecular docking using the GLIDE program (version 10.2, Schrodinger, LLC, New York, 2015). The crystal structure of the HDAC2 protein (PDB ID: 4LXZ; http://www.rcsb.org/pdb) was downloaded and the Protein Preparation Wizard module (Maestro, Schrodinger, LLC, New York, 2015) was applied for the structure preparation (i.e., align bond orders, added hydrogen, filled in missing side chains and loops using Prime, optimization with water molecules removed and then minimization with the OPLS_2005 force field, etc.). The native crystal structure of HDAC2 has three chains, A, B, and C, which all contain the same native ligand (Vorinostat) in different binding sites. Among the three protein chain-Vorinostat docking scores, the best score belongs to chain A. Therefore, chain A was selected for docking. The ligands were minimized by means of the LigPrep module (Schrödinger, LLC, New York, 2015). The candidate molecular ligands generated the proper ionization, tautomers, stereochemistries, and ring conformations. Partial

ΔG bind = Gcpx − (Grec + Glig)

(1)

The MM-PBSA (Molecular Mechanics-Poisson-Boltzmann/Surface Area) methodology allows the calculation of the complete binding reaction energy, including the desolvation of the ligand and the unbound protein, on the basis of a thermodynamic cycle. Therefore, Eq. (1) can be approximated as:

ΔG bind = ΔEMM –TΔS + ΔGsol

(2)

ΔEMM = ΔEele + ΔE vdw

(3)

All energies expressed in the above equations were averaged over the course of the molecular dynamics trajectories. In Eq. (3), ΔEMM is 65

Gene 637 (2017) 63–71

H. Zhou et al.

Fig. 2. Molecular structures of designed three series of urushiol derivatives (series 1: compounds 1–4; series 2: compounds 5–8; series 3: compounds 9–30).

nonpolar solvation term. The former component was calculated using the PB calculation, whereas the latter term is determined using Eq. (5):

the molecular mechanical energy obtained from the electrostatic (ΔEele) and the van der Waals (vdW, ΔEvdw) interactions within the system. Here, TΔS is the solute entropic contribution at temperature T (Kelvin) and the solvation free energy (ΔGsol) represents the electrostatic and nonpolar free energy of solvation, and therefore can be expressed as:

ΔGsol = ΔGele,sol + ΔGnp,sol

ΔGnp,sol = γ SASA + b

(5)

where SASA is the solvent-accessible surface area (Å ) and γ and b represent experimental solvation parameters. The normal mode calculation is required to obtain the entropy, but it is time-consuming for large systems, considering the entropic 2

(4)

where ΔGele,sol is the polar contribution to solvation and ΔGnp,sol is the 66

Gene 637 (2017) 63–71

H. Zhou et al.

2015). Wegener et al. reported that trifluoromethyl or trifluoromethoxy groups on benzene ring of TSA could increase its HDACI activity by two-fold (Wegener et al., 2008). Some studies have suggested that the introduction of amino, sulfonamido or acylamino groups at the cap moiety can significantly increase HDACI activity (Zhou et al., 2015). It has also been shown that the steric bulk of substituent groups could influence the HDACI activity; for instance, Zhang et al. reported phenylsulfonamido derivatives of SAHA demonstrated were more potent than its methylsulfonyl amino derivatives (Zhang et al., 2015). Furthermore, the introduction of electron-donating or -withdrawing or less bulky groups such as methyl, methoxyl, hydroxy or nitro at the cap moiety could also increase the activity of HDACIs (Nam et al., 2014). To further increase the binding affinity, selectivity and inhibition activity of urushiol, we introduced different substituent groups with different electron, steric bulk and hydrophobicity characteristics into the benzene ring of hydroxamic acid urushiol derivatives and explored the effects computationally.

differences among the systems is normally small when the receptors of all system models are quite similar, the entropy contribution of complexes was not considered in our ΔGbind calculations to save the computational cost; after all, our goal is to compare the binding energy of various ligands with the same receptor (as opposed to comparing the binding energies across receptors). 3. Results and discussion 3.1. Design rationale A cap group, a linker chain and a Zn2 + binding group are the essential pharmacophore for HDACIs (Yang et al., 2014). Molecular modeling studies showed that HDAC inhibitors appear to bind in a pocket of the HDAC active site in such a way that the metal binding moiety coordinates the catalytic ion at the end of the pocket, so the zinc binding group has proven to be crucial for the HDAC inhibition activity (Musso et al., 2015). The hydroxamic acid group is often used as a zinc binding group, which shows excellent chelation with Zn2 + by the two oxygen atoms in its structure (Chen et al., 2013). Urushiol has potent anticancer activities, and studies showed it also possess HDAC inhibition activity; however, urushiol has no zinc binding group, causing it to possess lower HDAC inhibition activity. Therefore, in our design, we first introduced the hydroxamic acid moiety into the structure of urushiol to make it possess the three portions of essential pharmacophore for HDACIs. Meanwhile, considering the ortho phenolic hydroxyl groups of urushiol can easily cause oxidative polymerization and the structural instability of urushiol, we carried out the chemical modification of the two adjacent phenolic hydroxyl groups by introducing different etherification or esterification groups. On the one hand, we expect to obtain structural stable urushiol derivatives, which will prevent the oxidative polymerization in the process of subsequent structural modification of urushiol. On the other hand, we expect to study the effects of the replacement of adjacent phenolic hydroxyl groups with different etherification or esterification groups on the HDAC inhibition activity of urushiol derivatives. The linker domain of the hydroxamic acid urushiol derivatives had a long alkyl chain with 10 carbon atoms. The linker domain of HDACIs plays an important role in their bioactivity (Li et al., 2012). There are numerous amino acid residues around the hydrophobic tunnel, and the introduction of some nitrogen or oxygen-containing groups into the linker domain can improve the binding effects through hydrogen bonding interactions (Liu et al., 2016). With these ideas in mind, we envisioned a strategy in which hydroxy, carbonyl, amino or methyl ether groups were introduced on the alkyl chain of the hydroxamic acid urushiol derivatives to occupy the HDAC binding pocket and confer better binding properties on the compounds. The surface recognition domain is another important composition for HDACIs, as it contributes significantly to the overall binding affinity and selectivity for HDAC (Pabba et al., 2011). The cap groups of the HDAC inhibitors can selectively interact with amino acid residues on the rim of the active pocket entrance of HDAC (Wang et al., 2014). Therefore, the cap groups of the inhibitors seem to be a promising region of design for improving the potency and selectivity of HDACIs. Many studies have demonstrated that the introduction of substituents in the cap group can significantly increase the binding affinity, selectivity, and specificity of HDACIs. Meanwhile, these substituents generating different electronic properties, steric hindrance or hydrophobic nature show different effects on the inhibitory activity toward HDAC (Clausen et al., 2015). Many studies show that nitrogen-containing heterocyclic groups (e.g., thiazole, pyrazole, pyrimidine, thiadiazole, triazole) placed at the surface recognition domain of HDACIs can significantly increase selectivity and inhibitory activity (Cai et al., 2015). It has been illustrated that the N, O and S atoms and the electron-rich nitrogencontaining heterocycle could form π–π and H-bonding interactions with the amino acid residues in the rim of the HDAC active site (Jin et al.,

3.2. Virtual screening The potency and selectivity of the designed 30 compounds were estimated by a grid-based scoring function to determine which of the urushiol derivatives could bind well to the active pocket of HDAC2 and thus have potential inhibitory activities. Studies have demonstrated that the Glide score is an excellent algorithm for screening active and inactive compounds. The Glide score is calculated from the contributions of hydrogen bonds and van der Waals between the protein and ligand, together with intra-molecular hydrogen bonds and ligand strains (Friesner et al., 2006). Table 1 shows the docking scores for the 30 urushiol derivatives. Compounds 21, 23, 10, 19, 9, 30, 1, 6, 5 and 29 gave the best scores of − 7.65 to − 8.47. These 10 compounds all fit the binding pocket quite well, the hydroxamic acid group formed a metal coordination with Zn2 +, and more than two hydrogen bonds formed with the residues. Hydrogen bond and hydrophobic interactions in the binding between the compounds and HDACs play an essential role in stabilizing the complex and increasing the binding affinity (Dixit et al., 2016). For the first series, compounds 1–4, compound 1 was found to display the best Glide score, and the methylene acetal group was found to be more potent as an electropositive group, possibly because there are chances of steric clashes if there is two phenylmethyl or methyl ethers or acetate ester groups. In the second series, structures 5–8, compounds 5 and 6 had the best Glide scores, suggesting that introducing an electron-donating group (hydroxy or carbonyl) into the alkyl side chain of urushiol is effective for increasing the binding affinity because they can form stronger hydrogen bond interactions with the residues. In the third series, structures 9–30, compounds 9, 10, 19, 21, 23, 29 and 30 exhibited the best Glide scores. These results indicated that addition of eF, eCl, triazole, benzamido, formamido, hydroxy or nitro substituents on the benzene ring could increase the docking scores significantly. This could be explained by that these groups could form π–π interaction and Table 1 Glide-scores of designed urushiol derivatives (compounds 1–30).

67

Compd.

Glide-score (kcal/mol− 1)

Compd.

Glide-score (kcal/mol− 1)

Compd.

Glide-score (kcal/mol− 1)

1 2 3 4 5 6 7 8 9 10

− 7.914 − 4.859 − 3.786 − 7.031 − 7.724 − 7.788 − 6.26 − 7.532 − 7.994 − 8.079

11 12 13 14 15 16 17 18 19 20

− 7.044 − 6.31 − 6.607 − 4.288 − 6.258 − 6.258 − 4.242 − 5.131 − 7.999 − 5.861

21 22 23 24 25 26 27 28 29 30

−8.474 −4.935 −8.178 −6.277 −5.123 −4.634 −7.627 −6.21 −7.655 −7.946

Gene 637 (2017) 63–71

H. Zhou et al.

Fig. 3. Molecular docking modeling of compounds in the active site of HDAC2 (a: HDAC2-compound 21; b: HDAC2compound 23; c: HDAC2-compound 10; d: HDAC2-compound 19; e: HDAC2-compound 9; f: HDAC2-compound 30).

phenyl group of compound 23 formed two π–π interactions with His183 and Phe155. Compound 19 formed a bidentate chelate with the active site Zn2 + through the triazole group on the benzene ring, and had hydrogen bonding interactions with Phe279, Gly273, Gly277, and Gly154 through the triazole and hydroxamate moieties. Moreover, the triazole and benzene ring of compound 19 formed five π–π interactions with His183, Tyr308, His146 and Phe155. It was found that when a triazole or amide moiety was introduced into the surface recognition domain, the carbonyl oxygen of amide group or the nitrogen of triazole group could form a chelate with Zn2 +, instead of the hydroxyl of hydroxamate group. Meanwhile, introducing a triazole or amide moiety could form more hydrogen bonding and π–π interactions. We hypothesize that the introduction of these groups into the cap structure of urushiol derivatives will enhance the interaction with the HDAC ligand binding pocket and further strengthen the HDACI activity.

H-bonding interaction with surrounding residues by the oxygen, nitrogen and halogen atoms in the groups. 3.3. Docking simulation study To gain some insights into the interaction between these compounds and HDAC2, six compounds (21, 23, 10, 19, 9, 30) with Glide scores better than SAHA were selected for docking simulation studies to elucidate the binding process, such as hydrogen bonding interactions and hydrophobic effects. As shown in Fig. 3, all compounds were successfully docked into the active pocket. Compounds 9, 10 and 30, showed a similar binding pattern with the aliphatic chain occupying the long and narrow channel, the hydroxamate group positioned at the bottom of this channel, the hydroxyl and carbonyl of hydroxamate in the three compounds could form hydrogen bonding interactions with His145 and Tyr308, and form a chelate with Zn2 +. Additionally, the imino group of compound 10 formed a hydrogen bond with His146, and the nitro of compound 30 formed hydrogen bonds with Glu103 and Asp104. The models for compounds 21 and 23, showed chelation of Zn2 + through the carbonyl oxygen of the benzamido or formamido moiety. Simultaneously, compound 21 had hydrogen bonding interactions with His183, Asp104 and Gly154 through its methylene acetal, hydroxamate and amide moiety, respectively. Compound 23 had hydrogen bonding interactions with Asp104, Glu103, and Gly154 through its hydroxamate and amide moiety. In addition to these hydrogen bonding interactions, the benzene methylene ether and benzamido moiety of compound 21 form four π–π interactions with His183, Phe155 and His146, and the

3.4. Molecular dynamics simulation To evaluate the stability and the flexibility of the docking complexes of HDAC2 with compounds 21, 23, 10, 19, 9 and 30 as a function of time, 20 ns MD simulations in an explicit hydration environment were performed. The results of MD simulations were analyzed in terms of stability of the trajectory, hydrogen bonds dynamics, and the conformational flexibility of HDAC2. Root mean square deviations (RMSD) of Cα, C, and N atoms for all complexes were analyzed to determine the stability of the trajectory during MD simulation. RMSD values 68

Gene 637 (2017) 63–71

H. Zhou et al.

Fig. 4. The RMSD of the compounds (21, 23, 10, 19, 9 and 30) during the 20 ns molecular dynamics.

Fig. 5. The RMSF for all residues during the last 500 ps in MD simulations.

residues in the active site, especially the residues at the bottom of the channel of HDAC2 (His145, His146, Gly154, Glu103, His183, Asp104, Tyr297, Gln254 and Tyr308) were all very small. This suggested that in HDAC2, some amino acids could shift far away from their normal positions when the enzyme is bound to an inhibitor, but for the amino acids that are in the active site, binding with an inhibitor could make them more rigid (Huang et al., 2012a). Hydrogen bonding between a ligand and a protein is very important for stabilizing the ligand-protein complex. In the current MD simulations, the hydrogen bonds between each of the compounds and HDAC2 were examined (Table 2). Such analysis showed that O and N atoms of compounds 21 and 23 were predicted to moderately hydrogen bond to Gly143 (occupancy 33.1 and 22.0%) and Glu92 (occupancy 14.3 and 68.4%); the O atoms of compounds 19 and 30 formed moderate hydrogen bonds to Ser342 and Asp93 with occupancy of 59.6 and 55.8%, respectively. Compound 10's O atoms strongly hydrogen bonded to Gly143, Tyr297 and His172, with occupancy of 65.1, 54.6, and 51.1%, respectively. Compound 9 formed weaker hydrogen bonds to His134 and Gln254, with occupancy of 18.4 and 29.1%. HDAC2 is a zinc-dependent enzyme; therefore, the ability of the compound to chelate the zinc ion is a precondition for the inhibition of the enzyme. The zinc ion coordination states in MD simulations showed that all the compounds could coordinate the zinc ion. In each complex, Zn2 + coordinated five or six atoms. The atoms included all carboxylate oxygens of Asp258 and Asp170, the nitrogen atom of His172, and the O

(compared to the energy-minimized starting structure) along the simulations are shown in Fig. 4. Small fluctuations and constant backbone atom RMSD are a good indication of system stability (Tanwar et al., 2014). As seen, after an initial period of RMSD fluctuations, all complexes reached equilibrium and remained stable during the remainder of the MD simulation, implying that the system folded to a more stable state than the starting structure. The complexes with compounds 21, 23 and 10 reached equilibrium after 5 ns of simulation, with average RMSD values of 1.78 Å, 1.41 Å, and 2.10 Å, respectively. The complexes with compounds 9, 19 and 30 reached equilibrium status after 1 ns, 4 ns and 8 ns, respectively, and their average RMSD values were 1.90 Å, 1.22 Å and 1.72 Å. The complex containing compound 9 showed a significant fluctuation within 15–16 ns, indicating conformational change due to spatial fitting of the ligand in the active site. After the fluctuation, it maintained a continuous equilibrium up to the end of the simulation time. These results suggest that the stabilities of the dynamic equilibriums for the complexes were reliable and that the trajectories could be useful in collecting snapshots for further analyses. In addition, the stability of the systems also bolstered the credibility of the docking results. Root mean squared fluctuation (RMSF) of backbone atoms for all complexes was analyzed to explore the flexibility of individual residues (Fig. 5). As shown in Fig. 5, all complexes had similar RMSF. Though some amino acid residues in some complexes fluctuated much more than in other complexes, the conformational changes of the amino acid 69

Gene 637 (2017) 63–71

H. Zhou et al.

Table 2 Hydrogen bonds between the compounds (21, 23, 10, 19, 9 and 30) and HDAC2 active site obtained from MD simulations. Compd.

Donor

Acceptor

Occupied(%)

Compd.

Donor

Acceptor

Occupied(%)

21

Compd.21:N2 Compd.21:O5 Compd.21:O5 Compd.23:N1 Compd.23:N2 Compd.23:O5 SER342:OG Compd.19:O4 GLY262:N Compd.19:N2

GLY143:O GLU92:OE2 GLN20:O GLU92:O GLY143:O GLN20:OE1 Compd.19:O3 ASN269:OD1 Compd.19:O4 GLY143:O

33.1 14.3 9.4 68.4 22.0 9.8 59.6 37.1 30.3 10.6

9

Compd.9:N1 Compd.9:O4 Compd.9:O4 Compd.9:O4 Compd.10:O4 TYR297:OH HIE172:NE2 Compd.30:O4 HIE172:NE2

HIE134:ND1 GLN254:OE1 GLY293:O ASP168:OD1 GLY143:O Compd.10:O3 Compd.10:O3 ASP93:OD2 Compd.30:O6

18.4 29.1 10.6 10.1 65.1 54.6 51.1 55.8 10.6

23

19

10

30

Table 3 Average distances between zinc ion and surrounding atoms in HDAC2-compound complexes. Compd.

Distance

21 23 10 19 9 30

Compd.atom

Asp258:OD1

Asp258:OD2

Asp170:OD1

Asp170:OD2

Hie172:ND1

Gly295:O

Compd.

1.69 1.85 1.78 1.78 1.78 2.07

1.86 1.75 1.82 1.83 1.83 1.73

1.79 1.88 1.88 1.86 1.87 1.78

1.84 1.78 1.79 1.76 1.77 1.85

2.13 2.06 2.12 – 2.11 –

– – – – – –

2.03 2.03 3.02 1.96 1.89 3.76

binding, whereas polar and non-polar desolvation energy terms (ΔGele,sol and ΔGnp,sol) opposed binding. In addition, by analyzing the non-polar interaction energies (ΔEvdw + ΔGnp,sol) and the polar interaction energies (ΔEele + ΔGele,sol), we found that the association between HDAC2 and compounds is mainly driven by a favorable nonpolar interaction (− 7.69 kcal mol− 1 to − 25.97 kcal mol− 1), whereas polar interactions are unfavorable to ligand binding (7.06 kcal mol− 1 to 27.61 kcal mol− 1). Meanwhile, the non-polar interaction energies of compounds 21, 23 and 19 (−25.97, − 20.04 and − 15.46 kcal mol− 1, respectively) were much higher than those of other compounds. It is hypothesized that the numerous aromatic groups of compounds 21, 23 and 19 and π–π stacking interactions with the aromatic residues of the receptor both contributed to this result.

or N atoms of the compounds, resulting in a pentahedral or hexahedral geometry (Table 3). The hydrogen bonds and Zn2 + coordination may be crucial for stabilizing the complexes, and this may be the reason for the potency of these compounds with respect to HDAC2 binding (Huang et al., 2012a).

3.5. Binding free energy calculations and energy decompositions To obtain more detailed information on the interactions between the compounds 21, 23, 10, 19, 9 and 30 and HDAC2, binding free energies were calculated by the MM/PBSA method. The results of the predicted binding free energies and energy components of each studied complex are listed in Table 4. Compounds 21 and 23 exhibited much lower binding energies of −6.88 and − 12.98 kcal mol− 1, respectively, compared to the binding energies of the rest of the compounds. Such negative values for binding energies indicate high binding affinity and thermodynamic stability (Sun et al., 2015). Compounds 10, 19, 9 and 30 showed binding energies of 2.35, 12.15, 2.91 and 2.45 kcal mol− 1, respectively. Such positive values of binding indicate less thermodynamic stability of the formed complexes. We can therefore predict that compounds 21 and 23 would have the best HDAC-inhibitory activity. This result was in accordance with their docking score results. To obtain a better view of which energy term had the greatest impact on the calculated binding affinities, the four individual energy components (ΔEvdw, ΔEele, ΔGele,sol and ΔGnp,sol) were compared. The data in Table 2 show that intermolecular van der Waals (ΔEvdw) and electrostatics (ΔEele) interactions provided major contributions to the

4. Conclusions In the present work, three series of 30 novel urushiol derivatives were designed as potential HDAC2 inhibitors. The design ideas included the following. A hydroxamic acid moiety was introduced into the tail of the alkyl side chain of urushiol and different ether or ester groups were introduced to replace the ortho phenolic hydroxyl groups. A hydroxy, carbonyl, amino or methyl ether group was introduced into the alkyl side chain of urushiol. Substituents of differing electronic or steric properties were introduced into the benzene ring. The binding affinity toward HDAC2 of all designed compounds was screened by Glide docking. Six compounds showed excellent Glide scores of − 7.9 to − 8.5, better than that of approved drug SAHA. The Glide docking studies revealed that introducing eF, eCl, triazole, benzamido, formamido, hydroxy or nitro substituents into the benzene ring could increase binding affinity significantly. Then, to further study the binding mode and mechanisms of the six best candidates as novel HDAC2 inhibitors, molecular docking, MD simulations, binding free energy calculations, and energy decomposition studies were carried out. Molecular docking studies revealed that Zn2 + coordination, hydrogen bonding, and hydrophobic interaction contributed to the high binding affinity of these compounds toward HDAC2. His145, His146, Gly154, Glu103, His183, Asp104, Tyr308 and Phe155 of HDAC2 contributed favorably to the binding between enzyme and the compounds. MD simulation and binding free energy calculations showed that all

Table 4 Binding free energies of compounds (21, 23, 10, 19, 9 and 30) to HDAC2 from MD simulations. Compd.

ΔEvdw

ΔEele

ΔGele,sol

ΔGnp,sol

ΔEMM

ΔGSOLV

ΔGbind

21 23 10 19 9 30

−56.29 −46.78 −32.52 −46.41 −38.46 −18.49

−56.74 −63.89 −36.48 −79.51 −64.60 −7.96

75.83 70.95 52.79 107.12 80.22 18.10

30.32 26.74 18.57 30.95 25.75 10.80

− 113.03 − 110.67 − 69.00 − 125.92 − 103.06 − 26.45

106.15 97.69 71.36 138.06 105.97 28.90

− 6.88 − 12.98 2.35 12.15 2.91 2.45

O1 O1 O3 N3 O3 O4

70

Gene 637 (2017) 63–71

H. Zhou et al.

cancer. Nat. Rev. Drug Discov. 1, 287–299. Li, X.H., Wei, Y.D., Wang, S.M., Wang, M.N., Huang, D.W., Xiu, Z.L., 2012. Synthesis, biological evaluation and molecular modeling of cyclic tetrapeptide based inhibitors HDAC. Chem. Res. Chin. Univ. 28, 1011–1016. Liu, R.S., Wang, J.H., Tang, W.P., Fang, H., 2016. Design and synthesis of a new generation of substituted purine hydroxamate analogs as histone deacetylase inhibitors. Bioorg. Med. Chem. 24, 1446–1454. Ma, X.M., Lu, R., Miyakoshi, T., 2012. Recent advances in research on lacquer allergy. Allergol. Int. 61, 45–50. Markiewitz, K.H., Dawson, C.R., 1965. On the isolation of the allegenically active components of the toxic principle of poison ivy. J. Org. Chem. 30, 1610–1613. Musso, L., Cincinelli, R., Zuco, V., Zunino, F., Nurisso, A., Cuendet, M., Giannini, G., Vesci, L., Pisano, C., Dallavalle, S., 2015. Investigation on the ZBG-functionality of phenyl-4-yl-acrylohydroxamic acid derivatives as histone deacetylase inhibitors. Bioorg. Med. Chem. Lett. 25, 4457–4460. Nam, N.H., Huong, T.L., Dung, D.T.M., Dung, P.T.P., Oanh, D.T.K., Park, S.H., Kim, K., Han, B.W., Yun, J., Kang, J.S., Kim, Y., Han, S.B., 2014. Synthesis, bioevaluation and docking study of 5-substitutedphenyl-1, 3, 4-thiadiazole-based hydroxamic acids as histone deacetylase inhibitors and antitumor agents. J. Enzyme Inhib. Med. Chem. 29, 611–618. Ning, C., Bi, Y., He, Y., Huang, W., Liu, L., Li, Y., Zhang, S., Liu, X., Yu, N., 2013. Design, synthesis and biological evaluation of di-substituted cinnamic hydroxamic acids bearing urea/thiourea unit as potent histone deacetylase inhibitors. Bioorg. Med. Chem. Lett. 23, 6432–6435. Pabba, C., Gregg, B.T., Kitchen, D.B., Chen, Z.J., Judkins, A., 2011. Design and synthesis of aryl ether and sulfone hydroxamic acids as potent histone deacetylase (HDAC) inhibitors. Bioorg. Med. Chem. Lett. 21, 324–328. Rikimaru, T., Taketomi, A., Yamashita, Y., Shirabe, K., Hamatsu, T., Shimada, M., Maehara, Y., 2007. Clinical significance of histone deacetylase 1 expression in patients with hepatocellular carcinoma. Oncology 72, 69–74. de Ruijter, A.J., van Gennip, A.H., Caron, H.N., Kemp, S., van Kuilenburg, A.B., 2003. Histone deacetylases (HDACs): characterization of the classical HDAC family. Biochem. J. 370, 737–749. Ryckewaert, L., Sacconnay, L., Carrupt, P.A., Nurisso, A., Simões-Pires, C., 2014. Nonspecific SIRT inhibition as a mechanism for the cytotoxicity of ginkgolic acids and urushiols. Toxicol. Lett. 229, 374–380. Sun, R., Li, X., Li, Y.Y., Zhang, X., Li, X.R., Li, X.Y., Shi, Z., Bao, J.K., 2015. Screening of novel inhibitors targeting lactate dehydrogenase a via four molecular docking strategies and dynamics simulations. J. Mol. Model. 21, 133. Tanwar, O., Deora, G.S., Tanwar, L., Kumar, G., Janardhan, S., Alam, M., Shaquiquzzaman, M., Akhter, M., 2014. Novel hydrazine derivatives as selective DPPIV inhibitors: findings from virtual screening and validation through molecular dynamics simulations. J. Mol. Model. 20, 2118. Wang, D., 2009. Computational studies on the histone deacetylases and the design of selective histone deacetylase inhibitors. Curr. Top. Med. Chem. 9, 241–256. Wang, B., Feig, M., Cukier, R.I., Burton, Z.F., 2013. Computational simulation strategies for analysis of multisubunit RNA polymerases. Chem. Rev. 113, 8546–8566. Wang, J.H., Sun, F.E., Han, L.Q., Hou, X.B., Pan, X.L., Liu, R.S., Tang, W.P., Fang, H., 2014. Design, synthesis, and preliminary bioactivity studies of substituted purine hydroxamic acid derivatives as novel histone deacetylase (HDAC) inhibitors. Med. Chem. Commun. 5, 1887–1891. Wang, C.Z., He, Y.F., Zhou, H., Tao, R., Chen, H.X., Ye, J.Z., Zhang, Y.S., 2015. Preparation and characterization of urushiol methylene acetal derivatives with various degrees of unsaturation in alkyl side chain. Int. J. Polym. Sci. 3, 1–6. Wegener, D., Deubzer, H.E., Oehme, I., Milde, T., Hildmann, C., Schwienhorst, A., Witt, O., 2008. HKI 46F08, a novel potent histone deacetylase inhibitor, exhibits antitumoral activity against embryonic childhood cancer cells. Anti-Cancer Drugs 19, 849–857. Xia, Z., Miyakoshi, T., Yoshida, T., 2004. Lipoxygenase-catalyzed polymerization of phenolic lipids suggests a new mechanism for allergic contact dermatitis induced by urushiol and its analogs. Biochem. Biophys. Res. Commun. 315, 704–709. Yamauchi, Y., Oshima, R., Kumanotani, J., 1980. Separation of Japanese lac urushiol diacetate on silver nitrate-coated silica gel columns by high-performance liquid chromatography. J. Chromatogr. A 198, 49–56. Yang, W., Li, L.X., Ji, X., Wu, X.W., Su, M.B., Sheng, L., Zang, Y., Li, J., Liu, H., 2014. Design, synthesis and biological evaluation of 4-anilinothieno [2, 3-d] pyrimidinebased hydroxamic acid derivatives as novel histone deacetylase inhibitors. Bioorg. Med. Chem. 22, 6146–6155. Yao, Y.W., Liao, C.Z., Li, Z., Wang, Z., Sun, Q., Liu, C.P., Yang, Y., Tu, Z.C., Jiang, S., 2014. Design, synthesis, and biological evaluation of 1, 3-disubstituted-pyrazole derivatives as new class I and IIb histone deacetylase inhibitors. Eur. J. Med. Chem. 86, 639–652. Zhang, S., Huang, W.B., Li, X.N., Yang, Z., Feng, B.H., 2015. Synthesis, biological evaluation, and computer-aided drug designing of new derivatives of hyperactive suberoylanilide hydroxamic acid histone deacetylase inhibitors. Chem. Boil. Drug. Des. 86, 795–804. Zhou, J.W., Li, M., Chen, N.H., Wang, S.L., Luo, H.B., Zhang, Y.K., Wu, R.B., 2015. Computational design of a time-dependent histone deacetylase 2 selective inhibitor. ACS Chem. Biol. 10, 687–692.

complexes possess good stability, which was characterized by low RMSD of the system, low RMSFs of residues, moderate hydrogen bonding, and Zn2 + coordination, and low values of binding free energies. In addition, the results of binding energy decomposition clearly showed van der Waals and electrostatics interactions provided major contributions to the stability of these complexes. This study shows the promising potential of urushiol derivatives as potent HDAC2 binding agents. Further studies on the synthesis and determination of their HDAC2 inhibitory properties are planned. Conflict of interest None. Acknowledgments This work was supported by the National Natural Science Foundation of China [grant number 31600467]; and the Jiangsu Key Lab. of Biomass Energy and Material Foundation of China [grant number JSBEM-S-201509]. References Abdel-Atty, M.M., Farag, N.A., Kassab, S.E., Serya, R.A., Abouzid, K.A., 2014. Design, synthesis, 3D pharmacophore, QSAR, and docking studies of carboxylic acid derivatives as histone deacetylase inhibitors and cytotoxic agents. Bioorg. Chem. 57, 65–82. Cai, J., Wei, H.T., Hong, K.H., Wu, X.Q., Zong, X., Cao, M., Wang, P., Li, L.S., Sun, C.L., Chen, B., Zhou, G.X., Chen, J.Q., Ji, M., 2015. Discovery, bioactivity and docking simulation of Vorinostat analogues containing 1, 2, 4-oxadiazole moiety as potent histone deacetylase inhibitors and antitumor agents. Bioorg. Med. Chem. 23, 3457–3471. Chen, G.L., Wang, L.H., Wang, J., Chen, K., Zhao, M., Sun, Z.Z., Wang, S., Zheng, H.L., Yang, J.Y., Wu, C.F., 2013. Discovery of a small molecular compound simultaneously targeting RXR and HADC: design, synthesis, molecular docking and bioassay. Bioorg. Med. Chem. Lett. 23, 3891–3895. Clausen, D.J., Smith, W.B., Haines, B.E., Wiest, O., Bradner, J.E., Williams, R.M., 2015. Modular synthesis and biological activity of pyridyl-based analogs of the potent class I histone deacetylase inhibitor largazole. Bioorg. Med. Chem. 23, 5061–5074. Dixit, V.A., Rathi, P.C., Bhagat, S., Gohlke, H., Petersen, R.K., Kristiansen, K., Chakraborti, A.K., Bharatam, P.V., 2016. Design and synthesis of novel Y-shaped barbituric acid derivatives as PPARγ activators. Eur. J. Med. Chem. 108, 423–435. Draper, W.M., Wijekoon, D., McKinney, M., Behniwal, P., Perera, S.K., Flessel, C.P., 2002. Atmospheric pressure ionization LC-MS-MS determination of urushiol congeners. J. Agric. Food Chem. 50, 1852–1858. Drummond, D.C., Noble, C.O., Kirpotin, D.B., Guo, Z., Scot, G.K.t, Benz, C.C., 2005. Clinical development of histone deacetylase inhibitors as anticancer agents. Annu. Rev. Pharmacol. Toxicol. 45, 495–528. Friesner, R.A., Murphy, R.B., Repasky, M.P., Frye, L.L., Greenwood, J.R., Halgren, T.A., Sanschagrin, P.C., Mainz, D.T., 2006. Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J. Med. Chem. 49, 6177–6196. Hess-Stumpp, H., 2005. Histone deacetylase inhibitors and cancer: from cell biology to the clinic. Eur. J. Cell Biol. 84, 109–121. Hildmann, C., Riester, D., Schwienhorst, A., 2007. Histone deacetylases-an important class of cellular regulators with a variety of functions. Appl. Microbiol. Biotechnol. 75, 487–497. Honda, T., Lu, R., Sakai, R., Ishimura, T., Miyakoshi, T., 2008. Characterization and comparison of Asian lacquer saps. Prog. Org. Coat. 61, 68–75. Hong, D.H., Han, S.B., Lee, C.W., Park, S.H., Jeon, Y.J., Kim, M.J., Kwak, S.S., Kim, H.M., 1999. Cytotoxicity of urushiols isolated from sap of Korean lacquer tree (Rhus vernicifera Stokes). Arch. Pharm. Res. 22, 638–641. Huang, D., Li, X., Sun, L., Xiu, Z., Nishino, N., 2012a. Synthesis, evaluation and molecular modeling of cyclic tetrapeptide histone deacetylase inhibitors as anticancer agents. J. Pept. Sci. 18, 242–251. Huang, D., Li, X., Wei, Y., Xiu, Z., 2012b. A novel series of l-2-benzyloxycarbonylamino-8(2-pyridyl)-disulfidyloctanoic acid derivatives as histone deacetylase inhibitors: design, synthesis and molecular modeling study. Eur. J. Med. Chem. 52, 111–122. Jin, K., Li, S.S., Li, X.G., Zhang, J., Xu, W.F., Li, X.C., 2015. Design, synthesis and preliminary biological evaluation of indoline-2, 3-dione derivatives as novel HDAC inhibitors. Bioorg. Med. Chem. 23, 4728–4736. Johnstone, R.W., 2002. Histone-deacetylase inhibitors: novel drugs for the treatment of

71