Original Research ARTICLE
Selective Inhibition of HDAC1 by Macrocyclic Polypeptide for the Treatment of Glioblastoma: A Binding Mechanistic Analysis Based on Molecular Dynamics
- 1College of Pharmaceutical Sciences, Zhejiang University, Hangzhou, China
- 2School of Pharmaceutical Sciences, Chongqing University, Chongqing, China
- 3State Key Laboratory of Applied Organic Chemistry and Department of Chemistry, Lanzhou University, Lanzhou, China
Glioblastoma (GBM) is the most common and aggressive intracranial malignant brain tumor, and the abnormal expression of HDAC1 is closely correlated to the progression, recurrence and metastasis of GBM cells, making selective inhibition of HDAC1 a promising strategy for GBM treatments. Among all available selective HDAC1 inhibitors, the macrocyclic peptides have gained great attention due to their remarkable inhibitory selectivity on HDAC1. However, the binding mechanism underlying this selectivity is still elusive, which increases the difficulty of designing and synthesizing the macrocyclic peptide-based anti-GBM drug. Herein, multiple computational approaches were employed to explore the binding behaviors of a typical macrocyclic peptide FK228 in both HDAC1 and HDAC6. Starting from the docking conformations of FK228 in the binding pockets of HDAC1&6, relatively long MD simulation (500 ns) shown that the hydrophobic interaction and hydrogen bonding of E91 and D92 in the Loop2 of HDAC1 with the Cap had a certain traction effect on FK228, and the sub-pocket formed by Loop1 and Loop2 in HDAC1 could better accommodate the Cap group, which had a positive effect on maintaining the active conformation of FK228. While the weakening of the interactions between FK228 and the residues in the Loop2 of HDAC6 during the MD simulation led to the large deflection of FK228 in the binding site, which also resulted in the decrease in the interactions between the Linker region of FK228 and the previously identified key amino acids (H134, F143, H174, and F203). Therefore, the residues located in Loop1 and Loop2 contributed in maintaining the active conformation of FK228, which would provide valuable hints for the discovery and design of novel macrocyclic polypeptide HDAC inhibitors.
Glioblastoma (GBM) is the most common and aggressive intracranial malignant brain tumor with the median survival duration <2 years in spite of chemotherapy, radiation or surgical resection (Natsume et al., 2019). In the current chemotherapies, such as temozolomide, drug resistance is the predominant obstacle (Zhang et al., 2012; Chen et al., 2019; Kim et al., 2019; Rahman et al., 2019; Su et al., 2019; Xingyi et al., 2019). On the basis of the latest experimental results obtained from the large-scale profiling which included the whole exome and RNA sequencing, it can be learnt that genetic and epigenetic mechanisms are involved in the occurrence and progress of glioma cells (Cancer Genome Atlas Research, 2008; Brennan et al., 2013), especially the aberrant epigenetic silencing of genes caused by histone deacetylation (Vaissiere et al., 2008; Cartron et al., 2013). A large number of researches have proven that significant nuclear expression of histone deacetylase 1 (HDAC1) occurred in GBM cells during the process of tumor progression, recurrence, and metastasis (Bhat et al., 2008; Kim et al., 2008; Campos et al., 2011; Li et al., 2016, 2018a; Zhang et al., 2016; Staberg et al., 2017; He et al., 2019; Natsume et al., 2019). In addition, the invasive and proliferative phenotype of GBM cells was found to be related to the overexpression of HDAC1 level (Han et al., 2013). Moreover, HDAC1 inhibitors developed for a variety of tumors have been extensively tested in clinical trials as a single drug or in combination with other chemotherapy agents (Lu et al., 2008; Tan et al., 2010; Campos et al., 2011; Dong et al., 2013; Tang et al., 2018). Currently, four HDAC inhibitors (HDACi) including Vorinostat, Romidepsin, Panobinostat, and Belinostat have been approved by FDA for anticancer therapeutics, and some other HDAC inhibitors (such as Ricolinostat) are still in the clinical trials to treat hematological and solid malignancies (Yang et al., 2016; Eckschlager et al., 2017; Li et al., 2018b).
Unfortunately, there are no clinical or approved cases of HDACi currently effective for the treatment of GBM. This is because targeting the key epigenetic enzymes, oncogenes, and pathways specific to glioblastoma cells by the drugs has proved to be of great challenges (Sturm et al., 2014), for example, the lower effective inhibitory concentrations within the tumor cells and adverse toxicological effects (Lee et al., 2015). In order to overcome the shortcoming caused by the limited stability and unacceptable pharmacokinetic properties of most existing drugs or molecules, various molecular skeletons were designed to improve the HDAC1-based drugs development, which conform the pharmacophore model of traditional HDACi, namely containing Cap group (Cap), Connect unit (CU), Linker region (Linker), and Zinc Binding Group (ZBG) (Figure 1; Dehmel et al., 2008; Varasi et al., 2011; Choi et al., 2012; Giannini et al., 2014; Krieger et al., 2017). Among these pharmacophores, the ZBG should penetrate deep into the bottom of the active pocket and chelate with zinc ion located in the catalytic center to compete with the protein for zinc ion, thereby inhibiting the catalytic activities of HDACs. And such binding pattern in the active pocket is the active conformation of the HDAC inhibitors (Krieger et al., 2019; Shen et al., 2019; Vergani et al., 2019).
Interestingly, the skeletons with macrocyclic Cap have better inhibitory activities against HDAC Class I than Class II, among which the macrocyclic peptide inhibitors account for a large proportion (Mwakwari et al., 2010; Rajak et al., 2013; Tapadar et al., 2015). As inhibiting class II HDACs (represented by HDAC6) can lead to unwanted toxic and side effects (especially serious cardiac toxicity) (Roche and Bertrand, 2016), targeting specific HDAC subtypes has shown great therapeutic potential. The macrocyclic HDACi targeting only Class I HDAC family or HDAC1 are regarded as lower toxicity and more tolerable than pan-HDAC inhibitors, which have shown great potential values of clinical therapeutic effects (Benelkebir et al., 2011; Bhansali et al., 2011; Mallinson and Collins, 2012; Salvador et al., 2014; Decroos et al., 2015; Pilon et al., 2015; Chen et al., 2017; Cheng et al., 2017; Kim et al., 2017). However, there are currently no crystal structures of HDAC1 and HDAC6 complexed with macrocyclic HDACi that have been resolved. Therefore, it is urgent and necessary to reveal the difference in the binding mechanism of macrocyclic HDACi in HDAC1&6 at the atomic level.
In this study, Romidepsin (FK228) was applied as a case study to investigate why macrocyclic polypeptide inhibitors tend to inhibit HDAC1, and various computational approaches were adopted to explore the binding modes of FK228 in HDAC1&6. First, the studied complexes of FK228 in HDAC1&6 were constructed via molecular docking approach. Second, the docked results were further verified by molecular dynamic (MD) simulation. Finally, the key residues responsible for the difference in binding energy of macrocyclic HDACi in HDAC1&6 were identified. In summary, the mechanism underlying why FK228 prefer to inhibit HDAC1 was elaborate through differential energy contributions and interaction fingerprints among the identified key amino acids, which could provide valuable information for the drug discovery on the basis of selective inhibition of HDAC1 in the future.
Results and Discussion
Construction of FK228 Complexed HDAC1and6 Structures
On the basis of the resolved protein crystals of HDAC1&6 available in Protein Data Bank (PDB) (Hai and Christianson, 2016; Watson et al., 2016), there were 75 and 68 binding poses of FK228 in HDAC1 and HDAC6 generated by molecular docking, respectively. Except for the docking score, the spatial similarity of the docking pose to the largazole thiol in HDAC8 was considered in selecting the conformations of FK228 in HDAC1&6 (Figure 2). This was because there were no HDAC1&6 crystals resolved with macrocyclic inhibitors, and HDAC8 (HDAC Class I) protein crystal with similar active pocket as HDAC1&6 complexed with largazole thiol could provide important clues for the choice of the initial conformations of FK228 in HDAC1&6 (Figure 2). According to Figure S1, it could be learnt that binding sites of HDAC1&6 were mainly composed of loop regions, namely loop 1–7. In addition, the selected docked poses suggested that the Cap group of FK228 had interactions with the residues at the rim of the active pocket of HDAC1&6, and the Linker coupled with ZBG penetrate the active pocket, which made the ZBG chelating with the zinc ion in catalytic center. Moreover, the orientation of FK228 in HDAC1&6 is highly coincident with the largazole thiol inhibitor in HDAC8 (Figure 2), which verified the reliability of the docking conformation to some extent. According to the Tables S1, S2, it could be found that the RMSD values were basically negatively correlated with the absolute value of the docking scores, and the smaller RMSD values could reflect the better binding of FK228 in the HDAC1&6 to some extent. In order to verify the reliability of the experiments, one additional initial conformation of the constructed system have been selected for the further molecular dynamic simulation (Figure S2).
Evaluating the Stability of MD Simulation via RMSD Analysis
The Complexes Stabilities Along the Simulation Monitored by RMSD
The selected docking conformations of HDAC1&6 in complex with FK228 were sampled by 500 ns MD simulation, and the dynamic trajectories of the studied complexes were supervised through the RMSD plots of the backbone-atoms of HDAC1&6, heavy-atoms of FK228, and the backbone-atoms of the amino acids in the binding pocket (within 5 Å of ligand) as the function of simulation time (Figure 3). Insight from the RMSD values in Figure 3, the FK228-HDAC1 and FK228-HDAC6 systems reached the equilibrium states around 350 and 50 ns, respectively. Moreover, according to the RMSD values of the additional independent simulations also showed that the constructed systems reached the equilibrium around 200 ns (Figure S3), and the difference in the fluctuation of the binding site of the two simulations were caused by the flexible loop domain.
Figure 3. Root mean square deviations of protein backbone atom, ligand heavy atoms, and the backbone atoms of the residues in the binding site as the function of time in MD simulations.
The Conformational Rearrangements of FK228 in HDAC1and6
The representative structures of FK228 binding to HDAC1&6 were obtained from the equilibrated trajectories and were compared with their corresponding initial conformations (Figure 4). During the MD process, the protein conformational changes were calculated by VMD software, the values were 1.85 and 1.92 Å for the HDAC1-FK228 and HDAC6-FK228 systems, indicating the small change in the conformation of protein. According to Figure 4A, it can be learnt that slight spatial shift of FK228 occurred in HDAC1 active site and the binding conformation maintain the interaction of sulfhydryl group (ZBG) chelating with the zinc ion (~3.2 Å) through inserting deeply into the active pocket. In contrast, for FK228 in HDAC6 (Figure 4B), there was a large deflection of the ZBG in the ligand from the initial conformation, namely straying from the catalytic center (~9.6 Å). In order to verify the reliability of the experiment, the conformational rearrangement of FK228 in HDAC1&6 of the additional independent experiment was also analyzed, and based on Figure S4, it could learnt that FK228 could maintain the active conformation in HDAC1 but not in HDAC6 (ZBG was also far away from the zinc ion). The conformational rearrangements investigated by MD simulation imply that the protein-ligand binding modes is the leading cause of the significant difference of FK228 inhibitory activity to HDAC1&6 and need to be further explored.
Figure 4. Comparison of the initial conformation and the representative conformation of the FK228 in HDAC1&6: (A) FK228 in HDAC1 system; (B) FK228 in HDAC6 system.
Molecular Mechanism of FK228 Selectivity to HDAC1and6
Insights From the FK228-HDAC1and6 Interaction Fingerprints
The binding modes of FK228 in HDAC1&6 are related to the interactions between drugs and amino acids of the target proteins. Thus, the interaction fingerprints analysis was used to explore the difference of FK228-HDAC1&6 binding modes (Figure 5). Figure 5A indicates that FK228 can maintain its interactions with the P22, E91, and D92 located at Loop1 and Loop2 of HDAC1 before and after MD simulation. For HDAC6-FK228 complex, although FK228 can maintain the interaction with P24 of Loop1, the interaction with S91 of Loop2 in the initial conformation was disappeared after MD simulation (Figure 5B). In addition, according to the interaction fingerprints, E91 locating on the Loop2 of HDAC1 contributed to a strong hydrophobic interaction with FK228, and the corresponding site on HDAC6 has no interaction with FK228, leading to the weak interaction between FK228 and Loop2 of HDAC6, which is the main reason of the large spatial shift of FK228 in the binding site of HDAC6.
Figure 5. Comparison of interaction fingerprints of FK228 in HDAC1&6 in the final 50 ns simulations with that of the optimized docking poses: (A) interaction fingerprints of FK228 in HDAC1; (B) interaction fingerprints of FK228 in HDAC6.
Insights From the Calculated Binding Free Energy of FK228-HDAC1and6 Complexes
The total binding free energies of HDAC1-FK228 and HDAC6-FK228 were −37.01 and −27.84 kcal/mol, which was consistent with the inhibitory gradient of FK228 toward HDAC1 and HDAC6 (Table 1). To qualify the energy contribution of each amino acid in HDAC1&6 for FK228's binding, the total binding free energies were decomposed at amino acid basis and the important ones with high contribution (≥0.1 kcal/mol) (Zheng et al., 2017) were identified. As shown in Figure 6 the values of amino acids energy with high contribution in each complex varied significantly (taking FK228-HDAC1 as example, the contribution of F143 equaled to −2.39 kcal/mol, which was almost 22 times of C93's energy contribution). As expected, the contributions of the amino acids at the corresponding position on HDAC1&6 also varied greatly. Taking G295 in HDAC1 and N306 in HDAC6 as example, it contributed −0.18 and −1.94 kcal/mol to the binding of FK228 in HDAC1 and HDAC6, respectively.
Table 1. Calculated and experimental data of FK228 binding to HDAC1 and HDAC6 (ΔG is in kcal/mol and IC50 value is in nM).
Figure 6. The per-residue binding free energy decomposition of 31 residues with high energy contribution (≥0.1 kcal/mol) to the interaction in at least one studied complex: FK228 in HDAC1 (light green); FK228 in HDAC6 (light orange).
In comparison with the residues D92 located at Loop2 of HDAC1, the reduction of the energy of corresponding residues S91 of HDAC6 contributed to FK228's binding enhanced our understanding of the difficulty of FK228 to maintain the initial conformation in HDAC6 of during MD simulation. As a result, the large spatial shift of FK228 in the binding site of HDAC6 led to the decreased energy contribution of the amino acids in the active site of HDAC6, such as F143, H174, and F203 (Figure 6), consisted very well with the decrease in the interacting frequency with these residues when compared with HDAC1 (Figure 7). For the Zn2+, the calculated energy contributing to FK228 binding in the active site of HDAC1 was −0.75 kcal/mol, while there was no energy contribution in HDAC6 system (Figure 6).
Figure 7. Comparison of the interaction fingerprints of FK228 in HDAC1&6 under the equilibrium trajectories.
The Key Role of Residue D92 in FK228 Binding to HDAC1
Both interaction fingerprints and amino acid energy contribution analysis found that residue D92 plays a key role in FK228 binding to HDAC1. The representative conformation obtained from the MD simulation trajectory in Figure 8 showed that the carbonyl group of D92 and nitrogen atom on the Cap group of FK228 could form a hydrogen bond. To further explore the huge difference in energy contribution of D92 in HDAC1 and its corresponding amino acid S91 in HDAC6, the distance between the two atoms forming the hydrogen bond during the equilibrium simulation (400–500 ns) was monitored, and the average distance between the two atoms forming hydrogen bonds was 3.15 Å (Figure 8). However, the side chain of S91 in HDAC6 lacked the hydrogen bond acceptor and its hydrophobic interaction with FK228 would gradually disappear with the deflection of its spatial position during the MD simulation (Figures 4B, 5B).
The Active Site Radius of Gyration Confirmed the Trend of FK228 to HDAC1
Physical and structural properties of the active pockets are closely related to the binding affinities of the ligands (Narang et al., 2019; Thillainayagam et al., 2019), the calculated binding free energy (Table 1) has successfully predicted the higher binding affinity of FK228 to HDAC1. To further evaluate the interactions between the FK228 and HDAC1&6, the radius of gyration (Rg) for the Ca-atoms of HDAC1&6 active pockets, which could be applied as an important and effective parameter to evaluate the structural integrity and compactness of the studied systems. The time evolution plot of Rg was calculated and shown in Figure 9. It is noted that the average Rg value of HDAC1-FK228 system is lower than that of HDAC6-FK228. The lower value of Rg of HDAC1-FK228 system indicated that the binding pocket of HDAC1 much more compacted and that FK228 could stay stably at the active site, which provided a guarantee for stronger interaction between FK228 and amino acids in the active pocket of HDAC1.
Overall Comparison of the Binding Conformations of FK228 in HDAC1and6
According to previous studies, HDAC inhibitors could exert the inhibitory activities by chelating with zinc ion at the catalytic center via the ZBG group deep into the bottom of the active pocket. According to Figure 10, the distance between the zinc ion and the sulfur atom on the ZBG of FK228 varied greatly in the two studied systems. In the HDAC1-FK228 system, the distance between the zinc ion and the sulfur atom on the ZBG of FK228 was about 3.5 Å, but the relative positions of sulfur and zinc ion is relatively larger in HDAC6-FK228 system. Furthermore, it can be learnt that Loop1 and Loop2 of HDAC1 formed a sub-pocket during the MD simulation process that could well-accommodate the sulfhydryl group on the Cap group and anchored the Cap group (Figure 11). The anchoring effects of Loop1 and Loop2 played a vital role in maintaining the binding conformation of FK228, and the relatively small spatial biases ensured the interaction of FK228 with important amino acids in the HDAC1 active pocket. However, for the HDAC6-FK228 system, the ZBG group did not penetrate the active pocket bottom to compete with the protein for metal zinc ions, and the previously calculated Rg value also indicated that the HDAC6 active pocket is less compacted, reducing the potential for interaction with FK228, which was the main reason for the large deflection of the Cap group at the active pocket.
Figure 10. Distance between the sulfur atom in the ZBG of FK228 and zinc ion in the studied systems.
As the first approved macrocyclic HDAC inhibitor, FK228 was used as a molecular probe to compare its binding conformation in HDAC1 and HDAC6, and to explore the molecular mechanism of FK228' tendency to inhibit HDAC1 at the atomic level through a variety of in silico approaches. For HDAC6-FK228 system, the disappearing hydrophobic interaction of S91 (located in the HDAC6 Loop2 region) with FK228 during the MD simulation and the lack of the corresponding residue of E91 (located in the HDAC1 Loop2 region) together weakened the anchoring effects of HDAC6 Loop2 to the FK228 Cap during the MD process, leading to the large spatial conformational deviation of the docking conformation and resulting the reduced interaction between the FK228 Linker region and the conserved amino acids in the HDAC6 pocket. In the case of HDAC1-FK228 system, K228 could maintain the interactions with D92 (hydrogen bonding) and E91 (hydrophobic interaction) on Loop2 after the dynamic trajectory reached equilibrium, and the interactions with H21 and P22 on Loop1 were also strengthened. Moreover, in the process of molecular dynamics simulation, Loop1 and Loop2 on HDAC1 could form a sub-pocket that better accommodated the Cap group of FK228, maintaining the active conformation of FK228 at the binding pocket and ensuring ZBG chelating with the zinc ion and competing with the protein for the metal zinc ion, thereby exerting the inhibitory activity on HDAC1. The interaction of the Cap group with the Loop1 and Loop2 regions contributes to maintaining the active conformation of the HDACi and should be especially considered on subsequent drug design based on selective inhibition of HDAC1.
Materials and Methods
The Construction of the Studied Systems
The studied systems FK228-HDAC1&6 were obtained by molecular docking using the Glide (2009) software embedded in Maestro (2009) with default parameters of standard precision. The 3D structure of FK228 was drawn by ChemBioDraw (Dickson et al., 2014) and saved in SDfile (*.sdf), then processed with LigPrep [OPLS-2005 (Price and Brooks, 2005) force fields] to generate the low-energy stable conformation. Additionally, the 3D structure of FK228 was preprocessed by Epik (2009) (pH = 7.0 ± 2.0) to generate the ionized state. After that, the protein structures of HDAC1&6 available in Protein Data Bank [PDB entry: 5ICN (Watson et al., 2016) and 5EDU (Hai and Christianson, 2016)] were processed by Protein Preparation Wizard (Maestro, 2009) module in Maestro (2009) to add the hydrogen atoms, assign protonation states and partial charges by OPLS-2005 (Price and Brooks, 2005) force field, and minimize the whole protein crystal to prepare the receptor for molecular docking. The minimization process is completed when the RMSD value reached 0.30 Å. Furthermore, the spatial coordinates of largazole analog in HDAC8 were referred when defining the docking grid due to the similar binding pockets (Cole et al., 2011; Du et al., 2011; Decroos et al., 2015; Gantt et al., 2016). In molecular docking, 5,000 poses were generated during the initial phase of the docking calculation, out of which best 400 poses were chosen for energy minimization by 100 steps of conjugate gradient minimizations (the details shown in Supplementary Materials).
Molecular Dynamics (MD) Simulation
MD simulation was performed within AMBER16 (2016) using GPU-accelerated PMEMD on 16 cores of an array of two 2.6 GHz Intel Xeon E5-2650v2 processors and 4 pieces of NVIDIA Tesla K40C graphic card. AMBER force field ff14SB (Dickson et al., 2014) and Li/Merz ion parameters (Li and Merz, 2014; Li et al., 2015a,b) were used for the protein and SPC/E water. General AMBER force field 2 (gaff2) was applied to assign the parameters of FK228 in each complex, and the atom types and partial charges of FK228 could be derived on the basis of RESP calculation through antechamber (Wang et al., 2006). In addition, the geometrical optimization and electrostatic potential calculation for FK228 were conducted at HF/6-31G* level through Gaussian 09 software (Gaussian 09, 2009). The Zn2+ was processed by 12-6-4 model (Li et al., 2015a) imbedded in Amber16. When the constructed systems were processed by LEaP (AMBER16, 2016), it could be found that FK228-HDAC1 and FK228-HDAC6 systems were solvated with a cubic water box, and the vdw box sizes were 527706.39 Å3 (12,126 water molecules) and 510534.33 Å3 (11,910 water molecules), respectively. In addition, there were two sodium ions in HDAC1-FK228 system that were used to neutralize the negative charge, and six sodium ions were used to neutralize the negative grid charge in the HDAC6-FK228 system.
Before the MD simulation, the processed research systems were subjected to the initial energy minimization through two procedure (Xue et al., 2018a; Zhang et al., 2019). The first step was to apply harmonic restraint on solute atom (force constant = 10 kcal·mol−1·Å−2), and the second step was to release all atoms to move freely. In each step, energy minimization was conducted by the steepest descent method for the first 5,000 steps and the conjugated gradient method for the subsequent 5,000 steps. Then, each studied system was heated from 0 to 100 K and then gradually to 310 K with the protein restrained over 100 ps in the NVT ensembles. Subsequently, 10 times (5 ns) unrestrained equilibration at 310 K were performed to equilibrate system's periodic boundary condition. Finally, the unrestrained 500 ns production simulation was conducted for the prepared four systems in NPT ensembles under the temperature of 310 K and the pressure of 1 atm. Temperature was controlled by Langevin dynamics and the pressure was controlled using Monte Carlo barostat (2016). In all the simulations, Particle-mesh Ewald (PME) (Darden et al., 1993) was used to handle the long range electrostatic interaction, and SHAKE algorithm was exploited to keep all bonds rigid (Larini et al., 2007; Fu et al., 2018; Xue et al., 2018a; Yang et al., 2018). Time step of simulation was set 2.0 fs and a 10.0 Å cutoff was used for non-bonded interactions (Xue et al., 2018b; Du et al., 2019).
All the analysis of MD trajectories, including as root mean square deviation (RMSD), the representative structures from the trajectories, binding free energies, were analyzed and predicted via cpptraj and mm_pbsa.pl programs embedded in AMBER16. Structural visualization was performed in PyMOL software (PyMOL 1.31).
Protein-Ligand Interaction Fingerprints Analysis
Interaction fingerprints between the FK228 and the HDAC1&6 were calculated via Ichem (Da Silva et al., 2018; Southan, 2018), and the calculation system mainly consisted of the ligand and the binding site (residues within 6 Å of the FK228's mass center). Firstly, conformation optimization and energy optimization were carried out for the docking poses of FK228 in HDAC1&6, and then interaction fingerprints was applied to carry out for the optimized conformations to calculate the interaction between FK228 and receptors in the initial conformation. Secondly, 500 snapshots were extracted from the equalized simulated trajectories (between 400 and 500 ns) to indicate the interacting effects between the FK228 and HDAC1&6, which was compared with interactions of the initial states of the studied systems. During the process of calculation, seven important interactions (hydrophobic interaction, aromatic, H-bond donor, H-bond acceptor, positively ionizable, negatively ionizable, and metal coordination) were applied to assess the interaction fingerprints between the ligand and receptor by parsing atoms and bond connectivity fields in the form of one-dimensional (1D) descriptors consisting of 1 and 0, and the results were shown by radar plots. In addition, detailed information about the rules of detecting the interactions between protein and ligand were shown in Table S3.
Calculation of the Binding Free Energy
MM/GBSA approach using a single molecular dynamic trajectory was adopted to calculate the binding free energy (ΔGMM/GBSA) regardless of entropic influence between the docked ligands and the receptor (Chen et al., 2017; Wang et al., 2017a, 2019; He et al., 2018), and in this study, 500 snapshots were extracted from the equilibrium trajectories (450–500 ns) for calculation via cpptraj. The calculation equation was as follows:
Where, ΔEvdW represented the van der Waals interactions contribution, ΔEele stood for the electrostatic energy contribution, ΔGpol was the polar solvent interaction energy calculated with the GB model (igb = 2) and Gnonpol was the non-polar solvation free energy, which was evaluated using LCPO method (0.0072 × ΔSASA, SASA indicating the solvent accessible area with a probe radius of 1.4 Å) (Weiser et al., 1999; Zheng et al., 2016).
Calculating the Per-Residue Energy Contribution
The per-residue energy contribution between the residues located in HDAC1&6 and the docked ligands was calculated using the following formula:
Radius of Gyration Calculation
The residues consisting of the binding site were selected to calculate the radius of gyration of the studied systems, and there were 30 residues in each systems. In this study, the equilibrium trajectories (450–500 ns) were used to calculate the Rg of the specified residues via cpptraj.
Data Availability Statement
All datasets generated for this study are included in the article/Supplementary Material.
FZ and WX conceived the work and directed the experiments. YZ performed the MD simulations. YZ, TF, YR, FL, GZ, JH, and XY collected and confirmed the data of protein and ligand structures and performed the analysis. FZ, WX, and YZ drafted the first and second version of the manuscript. All authors read, edited, and approved the final version of manuscript.
This study was supported by National Natural Science Foundation of China (81872798), National Key Research and Development Program of China (2018YFC0910500), Innovation Project on Industrial Generic Key Technologies of Chongqing (cstc2015zdcy-ztzx120003), and Fundamental Research Funds for the Central Universities (2018QNA7023, 10611CDJXZ238826, 2018CDQYSG0007, CDJZR14468801, and CDJKXB14011).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2020.00041/full#supplementary-material
Cap, capping group; CD, catalytic domain; CU, connect unit; GBM, Glioblastoma; HATs, histone acetyltransferases; HDAC6, Histone deacetylase 6; HDACi, HDAC inhibitors; MD, Molecular Dynamics; MM/GBSA, Molecular Mechanics Generalized Born Surface Area; sHDAC6Is, selectivity of HDAC6 inhibitors; ZBG, zinc binding group
1. ^PyMOL Molecular Graphics System, v. 1.3. Schrödinger, LLC, New York, NY.
Benelkebir, H., Marie, S., Hayden, A. L., Lyle, J., Loadman, P. M., Crabb, S. J., et al. (2011). Total synthesis of largazole and analogues: HDAC inhibition, antiproliferative activity and metabolic stability. Bioorg. Med. Chem. 19, 3650–3658. doi: 10.1016/j.bmc.2011.02.024
Bhansali, P., Hanigan, C. L., Casero, R. A., and Tillekeratne, L. M. (2011). Largazole and analogues with modified metal-binding motifs targeting histone deacetylases: synthesis and biological evaluation. J. Med. Chem. 54, 7453–7463. doi: 10.1021/jm200432a
Bhat, K. P., Pelloski, C. E., Zhang, Y., Kim, S. H., Delacruz, C., Rehli, M., et al. (2008). Selective repression of YKL-40 by NF-kappaB in glioma cell lines involves recruitment of histone deacetylase-1 and−2. FEBS Lett. 582, 3193–3200. doi: 10.1016/j.febslet.2008.08.010
Campos, B., Bermejo, J. L., Han, L., Felsberg, J., Ahmadi, R., Grabe, N., et al. (2011). Expression of nuclear receptor corepressors and class I histone deacetylases in astrocytic gliomas. Cancer Sci. 102, 387–392. doi: 10.1111/j.1349-7006.2010.01792.x
Cartron, P. F., Blanquart, C., Hervouet, E., Gregoire, M., and Vallette, F. M. (2013). HDAC1-mSin3a-NCOR1, Dnmt3b-HDAC1-Egr1 and Dnmt1-PCNA-UHRF1-G9a regulate the NY-ESO1 gene expression. Mol. Oncol. 7, 452–463. doi: 10.1016/j.molonc.2012.11.004
Chen, J. C., Lee, I. N., Huang, C., Wu, Y. P., Chung, C. Y., Lee, M. H., et al. (2019). Valproic acid-induced amphiregulin secretion confers resistance to temozolomide treatment in human glioma cells. BMC Cancer 19:756. doi: 10.1186/s12885-019-5843-6
Chen, X., Wang, C., Tang, S., Yu, C., and Zou, Q. (2017). CMSA: a heterogeneous CPU/GPU computing system for multiple similar RNA/DNA sequence alignment. BMC Bioinformatics 18:315. doi: 10.1186/s12859-017-1725-6
Cheng, K., Li, S., and Liao, C. (2017). Progress in the discovery of macrocyclic histone deacetylase inhibitors for the treatment of cancer. Curr. Med. Chem. 24, 4166–4179. doi: 10.2174/0929867324666170209105315
Choi, E., Lee, C., Cho, M., Seo, J. J., Yang, J. S., Oh, S. J., et al. (2012). Property-based optimization of hydroxamate-based gamma-lactam HDAC inhibitors to improve their metabolic stability and pharmacokinetic profiles. J. Med. Chem. 55, 10766–10770. doi: 10.1021/jm3009376
Cole, K. E., Dowling, D. P., Boone, M. A., Phillips, A. J., and Christianson, D. W. (2011). Structural basis of the antiproliferative activity of largazole, a depsipeptide inhibitor of the histone deacetylases. J. Am. Chem. Soc. 133, 12474–12477. doi: 10.1021/ja205972n
Decroos, C., Clausen, D. J., Haines, B. E., Wiest, O., Williams, R. M., and Christianson, D. W. (2015). Variable active site loop conformations accommodate the binding of macrocyclic largazole analogues to HDAC8. Biochemistry 54, 2126–2135. doi: 10.1021/acs.biochem.5b00010
Dehmel, F., Weinbrenner, S., Julius, H., Ciossek, T., Maier, T., Stengel, T., et al. (2008). Trithiocarbonates as a novel class of HDAC inhibitors: SAR studies, isoenzyme selectivity, and pharmacological profiles. J. Med. Chem. 51, 3985–4001. doi: 10.1021/jm800093c
Dong, L. H., Cheng, S., Zheng, Z., Wang, L., Shen, Y., Shen, Z. X., et al. (2013). Histone deacetylase inhibitor potentiated the ability of MTOR inhibitor to induce autophagic cell death in Burkitt leukemia/lymphoma. J. Hematol. Oncol. 6:53. doi: 10.1186/1756-8722-6-53
Du, J., Qiu, M., Guo, L., and Yao, X. (2019). Computational study of the binding mechanism between farnesoid X receptor alpha and antagonist N-benzyl-N-(3-(tertbutyl)-4-hydroxyphenyl)-2,6-dichloro-4-(dimethylamino) benzamide. J. Biomol. Struct. Dyn. 37, 1628–1640. doi: 10.1080/07391102.2018.1462735
Du, J., Sun, H., Xi, L., Li, J., Yang, Y., Liu, H., et al. (2011). Molecular modeling study of checkpoint kinase 1 inhibitors by multiple docking strategies and prime/MM-GBSA calculation. J. Comput. Chem. 32, 2800–2809. doi: 10.1002/jcc.21859
Fu, T., Zheng, G., Tu, G., Yang, F., Chen, Y., Yao, X., et al. (2018). Exploring the binding mechanism of metabotropic glutamate receptor 5 negative allosteric modulators in clinical trials by molecular dynamics simulations. ACS Chem. Neurosci. 9, 1492–1502. doi: 10.1021/acschemneuro.8b00059
Gantt, S. M., Decroos, C., Lee, M. S., Gullett, L. E., Bowman, C. M., Christianson, D. W., et al. (2016). General base-general acid catalysis in human histone deacetylase 8. Biochemistry 55, 820–832. doi: 10.1021/acs.biochem.5b01327
Giannini, G., Vesci, L., Battistuzzi, G., Vignola, D., Milazzo, F. M., Guglielmi, M. B., et al. (2014). ST7612AA1, a thioacetate-omega(gamma-lactam carboxamide) derivative selected from a novel generation of oral HDAC inhibitors. J. Med. Chem. 57, 8358–8377. doi: 10.1021/jm5008209
He, M., Li, W., Zheng, Q., and Zhang, H. (2018). A molecular dynamics investigation into the mechanisms of alectinib resistance of three ALK mutants. J. Cell. Biochem. 119, 5332–5342. doi: 10.1002/jcb.26666
Kim, B., Ratnayake, R., Lee, H., Shi, G., Zeller, S. L., Li, C., et al. (2017). Synthesis and biological evaluation of largazole zinc-binding group analogs. Bioorg. Med. Chem. 25, 3077–3086. doi: 10.1016/j.bmc.2017.03.071
Kim, B., Rincon Castro, L. M., Jawed, S., and Niles, L. P. (2008). Clinically relevant concentrations of valproic acid modulate melatonin MT(1) receptor, HDAC and MeCP2 mRNA expression in C6 glioma cells. Eur. J. Pharmacol. 589, 45–48. doi: 10.1016/j.ejphar.2008.04.058
Kim, H., Chong, K., Ryu, B. K., Park, K. J., Yu, M. O., Lee, J., et al. (2019). Repurposing penfluridol in combination with temozolomide for the treatment of glioblastoma. Cancers (Basel) 11:E1310. doi: 10.3390/cancers11091310
Krieger, V., Hamacher, A., Cao, F., Stenzel, K., Gertzen, C. G. W., Schaker-Hubner, L., et al. (2019). Synthesis of peptoid-based class I-selective histone deacetylase inhibitors with chemosensitizing properties. J. Med. Chem. 62, 11260–11279. doi: 10.1021/acs.jmedchem.9b01489
Krieger, V., Hamacher, A., Gertzen, C. G. W., Senger, J., Zwinderman, M. R. H., Marek, M., et al. (2017). Design, multicomponent synthesis, and anticancer activity of a focused histone deacetylase (HDAC) inhibitor library with peptoid-based cap groups. J. Med. Chem. 60, 5493–5506. doi: 10.1021/acs.jmedchem.7b00197
Larini, L., Mannella, R., and Leporini, D. (2007). Langevin stabilization of molecular-dynamics simulations of polymers by means of quasisymplectic algorithms. J. Chem. Phys. 126:104101. doi: 10.1063/1.2464095
Lee, P., Murphy, B., Miller, R., Menon, V., Banik, N. L., Giglio, P., et al. (2015). Mechanisms and clinical significance of histone deacetylase inhibitors: epigenetic glioblastoma therapy. Anticancer Res. 35, 615–625.
Li, P., Song, L. F., and Merz, K. M. Jr. (2015a). Parameterization of highly charged metal ions using the 12-6-4 LJ-type nonbonded model in explicit water. J. Phys. Chem. B 119, 883–895. doi: 10.1021/jp505875v
Li, S., Chen, X., Mao, L., Zahid, K. R., Wen, J., Zhang, L., et al. (2018a). Histone deacetylase 1 promotes glioblastoma cell proliferation and invasion via activation of PI3K/AKT and MEK/ERK signaling pathways. Brain Res. 1692, 154–162. doi: 10.1016/j.brainres.2018.05.023
Li, Y. H., Yu, C. Y., Li, X. X., Zhang, P., Tang, J., Yang, Q., et al. (2018b). Therapeutic target database update 2018: enriched resource for facilitating bench-to-clinic research of targeted therapeutics. Nucleic Acids Res. 46, D1121–D1127. doi: 10.1093/nar/gkx1076
Li, Z. Y., Li, Q. Z., Chen, L., Chen, B. D., Wang, B., Zhang, X. J., et al. (2016). Histone deacetylase inhibitor RGFP109 overcomes temozolomide resistance by blocking NF-kappaB-dependent transcription in glioblastoma cell lines. Neurochem. Res. 41, 3192–3205. doi: 10.1007/s11064-016-2043-5
Lu, Q., Lin, X., Feng, J., Zhao, X., Gallagher, R., Lee, M. Y., et al. (2008). Phenylhexyl isothiocyanate has dual function as histone deacetylase inhibitor and hypomethylating agent and can inhibit myeloma cell growth by targeting critical pathways. J. Hematol. Oncol. 9, 1–6. doi: 10.1186/1756-8722-1-6
Narang, S. S., Goyal, D., and Goyal, B. (2019). Inhibition of Alzheimer's amyloid-beta42 peptide aggregation by a bi-functional bis-tryptoline triazole: key insights from molecular dynamics simulations. J. Biomol. Struct. Dyn. 3, 1–14. doi: 10.1080/07391102.2019.1614093
Natsume, A., Hirano, M., Ranjit, M., Aoki, K., and Wakabayashi, T. (2019). Aberrant transcriptional regulation of super-enhancers by RET finger protein-histone deacetylase 1 complex in glioblastoma: chemoresistance to temozolomide. Neurol. Med. Chir. (Tokyo). 59, 293–298. doi: 10.2176/nmc.ra.2019-0049
Pilon, J. L., Clausen, D. J., Hansen, R. J., Lunghofer, P. J., Charles, B., Rose, B. J., et al. (2015). Comparative pharmacokinetic properties and antitumor activity of the marine HDACi Largazole and Largazole peptide isostere. Cancer Chemother. Pharmacol. 75, 671–682. doi: 10.1007/s00280-015-2675-1
Price, D. J., and Brooks, C. L. 3rd (2005). Detailed considerations for a balanced and broadly applicable force field: a study of substituted benzenes modeled with OPLS-AA. J. Comput. Chem. 26, 1529–1541. doi: 10.1002/jcc.20284
Rahman, M. A., Gras Navarro, A., Brekke, J., Engelsen, A., Bindesboll, C., Sarowar, S., et al. (2019). Bortezomib administered prior to temozolomide depletes MGMT, chemosensitizes glioblastoma with unmethylated MGMT promoter and prolongs animal survival. Br. J. Cancer 121, 545–555. doi: 10.1038/s41416-019-0551-1
Rajak, H., Singh, A., Dewangan, P. K., Patel, V., Jain, D. K., Tiwari, S. K., et al. (2013). Peptide based macrocycles: selective histone deacetylase inhibitors with antiproliferative activity. Curr. Med. Chem. 20, 1887–1903. doi: 10.2174/0929867311320140006
Salvador, L. A., Park, H., Al-Awadhi, F. H., Liu, Y., Kim, B., Zeller, S. L., et al. (2014). Modulation of activity profiles for largazole-based HDAC inhibitors through alteration of prodrug properties. ACS Med. Chem. Lett. 5, 905–910. doi: 10.1021/ml500170r
Shen, S., Hadley, M., Ustinova, K., Pavlicek, J., Knox, T., Noonepalle, S., et al. (2019). Discovery of a new isoxazole-3-hydroxamate-based histone deacetylase 6 inhibitor SS-208 with antitumor activity in syngeneic melanoma mouse models. J. Med. Chem. 62, 8557–8577. doi: 10.1021/acs.jmedchem.9b00946
Staberg, M., Michaelsen, S. R., Rasmussen, R. D., Villingshoj, M., Poulsen, H. S., and Hamerlik, P. (2017). Inhibition of histone deacetylases sensitizes glioblastoma cells to lomustine. Cell. Oncol. (Dordr). 40, 21–32. doi: 10.1007/s13402-016-0301-9
Sturm, D., Bender, S., Jones, D. T., Lichter, P., Grill, J., Becher, O., et al. (2014). Paediatric and adult glioblastoma: multiform (epi)genomic culprits emerge. Nat. Rev. Cancer 14, 92–107. doi: 10.1038/nrc3655
Tang, W., Wan, S., Yang, Z., Teschendorff, A. E., and Zou, Q. (2018). Tumor origin detection with tissue-specific miRNA and DNA methylation markers. Bioinformatics 34, 398–406. doi: 10.1093/bioinformatics/btx622
Tapadar, S., Fathi, S., Raji, I., Omesiete, W., Kornacki, J. R., Mwakwari, S. C., et al. (2015). A structure-activity relationship of non-peptide macrocyclic histone deacetylase inhibitors and their anti-proliferative and anti-inflammatory activities. Bioorg. Med. Chem. 23, 7543–7564. doi: 10.1016/j.bmc.2015.10.045
Thillainayagam, M., Ramaiah, S., and Anbarasu, A. (2019). Molecular docking and dynamics studies on novel benzene sulfonamide substituted pyrazole-pyrazoline analogues as potent inhibitors of Plasmodium falciparum Histo aspartic protease. J. Biomol. Struct. Dyn. 24, 1–11. doi: 10.1080/07391102.2019.1654923
Varasi, M., Thaler, F., Abate, A., Bigogno, C., Boggio, R., Carenzi, G., et al. (2011). Discovery, synthesis, and pharmacological evaluation of spiropiperidine hydroxamic acid based derivatives as structurally novel histone deacetylase (HDAC) inhibitors. J. Med. Chem. 54, 3051–3064. doi: 10.1021/jm200146u
Vergani, B., Sandrone, G., Marchini, M., Ripamonti, C., Cellupica, E., Galbiati, E., et al. (2019). Novel benzohydroxamate-based potent and selective histone deacetylase 6 (HDAC6) inhibitors bearing a pentaheterocyclic scaffold: design, synthesis, and biological evaluation. J. Med. Chem. 62, 10711–10739. doi: 10.1021/acs.jmedchem.9b01194
Wang, E., Sun, H., Wang, J., Wang, Z., Liu, H., Zhang, J. Z. H., et al. (2019). End-point binding free energy calculation with MM/PBSA and MM/GBSA: strategies and applications in drug design. Chem Rev. 119, 9478–9508. doi: 10.1021/acs.chemrev.9b00055
Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2006). Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 25, 247–260. doi: 10.1016/j.jmgm.2005.12.005
Wang, P., Fu, T., Zhang, X., Yang, F., Zheng, G., Xue, W., et al. (2017a). Differentiating physicochemical properties between NDRIs and sNRIs clinically important for the treatment of ADHD. Biochim. Biophys. Acta 1861, 2766–2777. doi: 10.1016/j.bbagen.2017.07.022
Wang, P., Zhang, X., Fu, T., Li, S., Li, B., Xue, W., et al. (2017b). Differentiating physicochemical properties between addictive and nonaddictive ADHD drugs revealed by molecular dynamics simulation studies. ACS Chem. Neurosci. 8, 1416–1428. doi: 10.1021/acschemneuro.7b00173
Watson, P. J., Millard, C. J., Riley, A. M., Robertson, N. S., Wright, L. C., Godage, H. Y., et al. (2016). Insights into the activation mechanism of class I HDAC complexes by inositol phosphates. Nat. Commun. 7:11262. doi: 10.1038/ncomms11262
Weiser, J., Peter, S., Shenkin, and Still, W. C. (1999). Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 20, 217–230. doi: 10.1002/(SICI)1096-987X(19990130)20:2217::AID-JCC43.0.CO;2-A
Xingyi, J., Guonan, C., Xin, Z., and Naijie, L. (2019). AbCD133 modified alphaCT1 loaded target magnetic mesoporous silica nano-drugcarriers can sensitizes glioma cancer stem cells to TMZ and have therapeutic potential on TMZ resistant glioblastoma. J. Biomed. Nanotechnol. 15, 1468–1481. doi: 10.1166/jbn.2019.2795
Xue, W., Wang, P., Tu, G., Yang, F., Zheng, G., Li, X., et al. (2018a). Computational identification of the binding mechanism of a triple reuptake inhibitor amitifadine for the treatment of major depressive disorder. Phys. Chem. Chem. Phys. 20, 6606–6616. doi: 10.1039/C7CP07869B
Xue, W., Yang, F., Wang, P., Zheng, G., Chen, Y., Yao, X., et al. (2018b). What contributes to serotonin-norepinephrine reuptake inhibitors' dual-targeting mechanism? the key role of transmembrane domain 6 in human serotonin and norepinephrine transporters revealed by molecular dynamics simulation. ACS Chem Neurosci. 9, 1128–1140. doi: 10.1021/acschemneuro.7b00490
Yang, F., Zheng, G., Fu, T., Li, X., Tu, G., Li, Y. H., et al. (2018). Prediction of the binding mode and resistance profile for a dual-target pyrrolyl diketo acid scaffold against HIV-1 integrase and reverse-transcriptase-associated ribonuclease H. Phys. Chem. Chem. Phys. 20, 23873–23884. doi: 10.1039/C8CP01843J
Yang, H., Qin, C., Li, Y. H., Tao, L., Zhou, J., Yu, C. Y., et al. (2016). Therapeutic target database update 2016: enriched resource for bench to clinical drug target and targeted pathway information. Nucleic Acids Res. 44, D1069–1074. doi: 10.1093/nar/gkv1230
Yurek-George, A., Cecil, A. R., Mo, A. H., Wen, S., Rogers, H., Habens, F., et al. (2007). The first biologically active synthetic analogues of FK228, the depsipeptide histone deacetylase inhibitor. J. Med. Chem. 50, 5720–5726. doi: 10.1021/jm0703800
Zhang, Y., Ying, J. B., Hong, J. J., Li, F. C., Fu, T. T., Yang, F. Y., et al. (2019). How does chirality determine the selective inhibition of histone deacetylase 6? A lesson from trichostatin a enantiomers based on molecular dynamics. ACS Chem. Neurosci. 10, 2467–2480. doi: 10.1021/acschemneuro.8b00729
Zhang, Z., Wang, Y., Chen, J., Tan, Q., Xie, C., Li, C., et al. (2016). Silencing of histone deacetylase 2 suppresses malignancy for proliferation, migration, and invasion of glioblastoma cells and enhances temozolomide sensitivity. Cancer Chemother. Pharmacol. 78, 1289–1296. doi: 10.1007/s00280-016-3188-2
Zheng, G., Xue, W., Wang, P., Yang, F., Li, B., Li, X., et al. (2016). Exploring the inhibitory mechanism of approved selective norepinephrine reuptake inhibitors and reboxetine enantiomers by molecular dynamics study. Sci. Rep. 6:26883. doi: 10.1038/srep26883
Zheng, G., Xue, W., Yang, F., Zhang, Y., Chen, Y., Yao, X., et al. (2017). Revealing vilazodone's binding mechanism underlying its partial agonism to the 5-HT1A receptor in the treatment of major depressive disorder. Phys. Chem. Chem. Phys. 19, 28885–28896. doi: 10.1039/C7CP05688E
Keywords: HDAC, macrocyclic peptides, molecular docking, MD simulation, binding free energies, interaction fingerprints
Citation: Zhang Y, Fu T, Ren Y, Li F, Zheng G, Hong J, Yao X, Xue W and Zhu F (2020) Selective Inhibition of HDAC1 by Macrocyclic Polypeptide for the Treatment of Glioblastoma: A Binding Mechanistic Analysis Based on Molecular Dynamics. Front. Mol. Biosci. 7:41. doi: 10.3389/fmolb.2020.00041
Received: 20 November 2019; Accepted: 21 February 2020;
Published: 11 March 2020.
Edited by:Hongchun Li, Shenzhen Institutes of Advanced Technology (CAS), China
Reviewed by:Quan Zou, University of Electronic Science and Technology of China, China
Sophie Sacquin-Mora, UPR9080 Laboratoire de Biochimie Théorique (LBT), France
Xingcheng Lin, Massachusetts Institute of Technology, United States
Copyright © 2020 Zhang, Fu, Ren, Li, Zheng, Hong, Yao, Xue and Zhu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.