Investigating the binding mechanism of AML inhibitors based on panobinostat with HDAC3 proteins using Gaussian accelerated molecular dynamics

Abstract

Class I histone deacetylases (HDACs) play a crucial role in the transformation and survival of myeloid and lymphoid malignancies, with HDAC1, 2, and 3 (Class I HDACs) being potential molecular targets for acute myelogenous leukemia (AML) treatment. Among them, HDAC3 depletion or inhibition significantly reduces proliferation and promotes differentiation in leukemia, with inhibitors like Panobinostat and compound 13a showing promise in suppressing its activity. In this study, we utilized Gaussian accelerated molecular dynamics (GaMD) simulations to compare the inhibitory potency of 13a and Panobinostat against HDAC3. Our findings suggest that the superior inhibitory activity of 13a may be attributed to its stronger interactions with HDAC3. Distance analysis demonstrated that 13a maintains a closer and more consistent coordination with the zinc ion in the catalytic pocket, resulting in a more stable interaction compared to Panobinostat. Additionally, interaction analysis revealed that 13a forms more π-alkyl interactions, along with additional attractive charge and metal–acceptor interactions with HDAC3. Principal component analysis (PCA) further showed that the binding of 13a stabilizes HDAC3 in multiple distinct conformational states, suggesting that a more substantial conformational rearrangement is required upon 13a binding. This structural complexity may explain why 13a behaves as a slow-on/slow-off inhibitor and exhibits a superior IC₅₀ compared to Panobinostat. Alanine scanning identified residues such as PRO23, HIS125, and PHE144 as potential sites for inhibitor binding, making significant contributions to binding affinity. These combined findings suggest that 13a not only has a higher inhibitory potency but also holds potential for further optimization, making it a promising candidate for targeted cancer therapy.

---

Introduction

Epigenetic dysregulation is a crucial factor in the onset and progression of cancer, as it can lead to the activation of oncogenes and the suppression of tumor suppressor genes.¹ These epigenetic changes are intricately involved in cancer development and metastasis.² A key aspect of epigenetic regulation is the post-translational acetylation of lysine residues on histone tails, a process mediated by histone acetyltransferases. This acetylation results in chromatin remodeling, which facilitates transcriptional activation.² On the other hand, the removal of these acetyl groups by histone deacetylases (HDACs) causes chromatin to become more compact, leading to changes in gene expression that can promote cancer progression.^3,4 As a result, HDACs have become attractive targets for the development of new anticancer therapies.

HDACs comprise a family of zinc-dependent metalloenzymes that are deeply involved in regulating cell migration and invasion in various types of cancer. The 18 isoforms of mammalian HDACs are categorized based on their structural, functional, and evolutionary characteristics into four classes: Class I (HDAC 1, 2, 3, and 8), Class II (HDAC 4, 5, 6, 7, 9, and 10), Class III (Sirtuins 1–7), and Class IV (HDAC 11).⁵ Among these, HDAC3, a member of Class I, is often overexpressed in a variety of cancers, such as acute myeloid leukemia (AML), prostate cancer, melanoma, breast cancer, and others.^6–10 As a key epigenetic regulator of multiple cellular signaling pathways, it plays a central role in cancer progression.^11,12

Acute myeloid leukemia (AML) is the most common form of adult leukemia and has the highest mortality rate among all types of leukemia in the United States.¹³ In the context of acute leukemia, HDAC3 activity is particularly important, as it is essential for the initiation of leukemia development and contributes to chemotherapy resistance by regulating DNA damage repair mechanisms.^14,15 Previous research has demonstrated that HDAC3 inhibitors can induce apoptosis in AML cells by indirectly inhibiting the FLT3/STAT5 signaling pathway and downregulating key anti-apoptotic proteins such as c-FLIP and XIAP.¹⁶ These findings underscore the potential of HDAC3 as a therapeutic target in AML and possibly other malignancies.

The FDA has approved four histone deacetylase inhibitors (HDACIs): Vorinostat,¹⁷ Romidepsin,¹⁸ Belinostat,¹⁹ and Panobinostat.²⁰ Among them, Panobinostat is the most potent HDACI in both in vitro and in vivo studies. The structural domains and pharmacophore model of HDACIs typically consist of three key components: a cap group (which interacts with the enzyme surface), a linker group (which occupies the long hydrophobic channel), and a zinc-binding group (ZBG, which acts at the catalytic site),^21,22 as shown in Fig. S1.† Hydroxamic acid is the most commonly used ZBG but has drawbacks such as structural instability, poor isomer selectivity, and mutagenicity. In recent years, hydrazide-based ZBG-HDACIs have emerged, aiming to explore new ZBGs that align with the geometry of physiological substrates as new scaffolds for HDACIs. However, their pharmacological effects, structural stability, and off-target toxicity still require further investigation.^23,24

C. James Chou and colleagues developed a series of HDAC inhibitors derived from Panobinostat's structure.²⁵ The lead compound, 13a, exhibited potent and selective inhibition of HDAC3 with an IC₅₀ of 0.28 nM, and displayed a unique noncompetitive binding mode. In wt-p53 FLT3-ITD MV4-11 leukemia cells, 13a downregulated FLT3, STAT5, and pERK, suppressed anti-apoptotic proteins such as XIAP and c-FLIP, and activated pro-caspase3, inducing apoptosis rather than autophagy. Its activity was p53- and FLT3-dependent, showing limited efficacy in p53-null or FLT3-wild-type cells. Compared to Panobinostat, 13a was non-mutagenic, did not induce Hsp70, and demonstrated improved pharmacokinetics including enhanced bioavailability, longer half-life, and a slow-on/slow-off inhibition profile, highlighting its potential as a selective and durable therapeutic agent for AML.²⁵ Despite these promising results, the precise inhibitory mechanism of 13a remains to be elucidated.

Gaussian accelerated molecular dynamics (GaMD) is an enhanced sampling method that accelerates conformational transitions by adding a harmonic boost potential to smooth the energy landscape. Considering the widespread application of GaMD in studying cancer-related protein–inhibitor interactions,^26–29 and given that 13a is a slow-on/slow-off inhibitor, we employed this enhanced sampling technique to investigate its binding mechanism with HDAC3. Specifically, we conducted simulations on three systems: HDAC3 alone, HDAC3 in complex with Panobinostat (HDAC3–Pan), and HDAC3 in complex with compound 13a (HDAC3–13a). Our goal was to elucidate the conformational changes induced in HDAC3 upon ligand binding and to explore the inhibitory mechanisms of di-N-substituted hydrazide-based HDACIs. The findings from this study may provide valuable insights into the dynamics of HDAC3 and contribute to the rational design of drugs for cancer therapy, ultimately improving therapeutic outcomes.

Results

Molecular docking analysis and system setup

We obtained the HDAC3 protein structure (PDB ID: 4A69) from the Protein Data Bank³⁰ and performed molecular docking with two inhibitors using AutoDock Vina.³¹ Similar to the study by C. James Chou et al.,²⁵ we retained the nuclear receptor corepressor 2 (NCoR2) chain from the 4A69 structure in the docking process. The 2D structures of Panobinostat (Pan) and 13a are shown in Fig. S2,† and the docking results are presented in Fig. 1. As shown in Fig. 1, both Pan and 13a exhibit similar binding modes, consistent with the typical binding pattern of HDAC inhibitors: the cap group interacts with the enzyme surface, the linker occupies the hydrophobic channel, and the zinc-binding group coordinates at the catalytic site.

Fig. 1 (A) The results of molecular docking correspond to system HDAC3–Pan. (B) The results of molecular docking correspond to system HDAC3–13a.

Docking results revealed that the binding energy of HDAC3–Pan was −10.99 kcal mol⁻¹, while that of HDAC3–13a was −14.87 kcal mol⁻¹, indicating that 13a has a stronger binding affinity for HDAC3 than Pan. Fig. S3 and S4† illustrates the interactions of HDAC3–Pan and HDAC3–13a. Pan forms hydrogen bonds with GLY21, ASP93, HIS134, HIS135, and ASP170, pi–pi interactions with PHE144 and PHE200, pi–cation interaction with HIS22, and pi–alkyl interactions with PRO23 and LEU266. Additionally, it engages in metal–acceptor interactions with zinc ion and van der Waals interactions with nearby residues such as ASP92, GLY143, and GLY296. 13a forms hydrogen bonds with ASP92, ASP93, HIS135, pi–pi interaction with PHE144, pi–cation interaction with HIS22, and pi–alkyl interactions with LEU266. Additionally, it engages in van der Waals interactions with nearby residues such as PRO23, GLY143, and GLY296. The stronger binding affinity of 13a may be attributed to the shorter and more stable hydrogen bonds it forms with HDAC3. Furthermore, 13a forms van der Waals interactions with more residues compared to Pan.

After determining the docking positions as the initial conformations, a 50 ns conventional MD (cMD) simulation was performed to pre-equilibrate the systems. As shown in Fig. S5,† after the initial equilibration period, all three systems—HDAC3, HDAC3–Pan, and HDAC3–13a—gradually stabilize in their respective conformations. The root mean square deviation (RMSD) values initially show greater variation but gradually converge to a more consistent and lower fluctuation pattern over time. This trend suggests that the proteins achieved a more stable structure, likely reflecting successful adaptation of the protein–ligand complexes to the simulation environment. Based on these results, further GaMD simulations were conducted. To validate the docking results, the representative structures obtained from k-means clustering of the GaMD trajectories were subsequently used for molecular docking with the inhibitors Pan and 13a. The binding energy results showed that HDAC3–Pan had a binding energy of −9.95 kcal mol⁻¹, whereas HDAC3–13a exhibited a lower binding energy of −12.34 kcal mol⁻¹. Consistent with the initial docking results, 13a maintained a stronger binding affinity.

Dynamic properties of the three systems

RMSD, which measures the average distance between protein atoms and their reference state after superposition, is widely used to assess the stability and conformational changes of molecular structures during simulations.^32,33 We calculated the RMSD for all three systems, as illustrated in Fig. 2A.

Fig. 2 (A) The RMSD diagrams of HADC3 and HDAC3–Pan, HDAC3–13a. (B) The R_g diagrams of HADC3 and HDAC3–Pan, HDAC3–13a. (C) The SASA diagrams of HADC3 and HDAC3–Pan, HDAC3–13a.

For the apo HDAC3 system (gray line), representing the unbound form of HDAC3, the RMSD gradually increases during the first 250 ns, starting at approximately 1 Å and reaching around 2 Å. Subsequently, it fluctuates between 1.5 and 2 Å.

In contrast, the HDAC3–Pan system (orange line) shows a rapid increase in RMSD within the first 50 ns. After this initial rise, the RMSD fluctuates between 1 and 2 Å, maintaining lower and more consistent values throughout the simulation. This indicates that the complex exhibits greater conformational stability compared to the apo form.

The HDAC3–13a system (blue line) displays the highest and most variable RMSD values, pointing to significant conformational changes throughout the simulation. The RMSD peaks multiple times, suggesting that the molecule frequently adopts different conformations. This pattern may indicate a more flexible structure or instability in its interaction with the 13a ligand.

Radius of gyration (R_g) serves as an important measure of a biomolecule's folding state, with lower values indicating greater compactness and higher values suggesting increased unfolding or irregularity in structure. By calculating R_g, we can observe the changes in compactness of the three systems during the simulation process.

As shown in Fig. 2B, both HDAC3–Pan and HDAC3–13a exhibit noticeably lower R_g values compared to the apo HDAC3 system, indicating that inhibitor binding induces a more compact protein conformation. These results suggest that 13a acts similarly to Pan, with both compounds promoting structural tightening of HDAC3. However, HDAC3–13a displays an even lower R_g than HDAC3–Pan, suggesting a greater degree of compaction.

The solvent-accessible surface area (SASA) is a key parameter used to evaluate the extent to which a molecule's surface is exposed to the surrounding solvent. It reflects the protein's interaction with its external environment and is sensitive to dynamic structural changes. Analyzing SASA changes during the simulation helps assess conformational shifts and surface exposure across the three systems.

As shown in Fig. 2C, the SASA values for the HDAC3 system remain relatively high throughout the 500 ns simulation, indicating greater surface exposure to solvent. In contrast, the SASA values of HDAC3–Pan and HDAC3–13a are noticeably lower than those of HDAC3, with HDAC3–13a showing the most substantial reduction throughout the simulation.

Flexibility and conformational dynamics analysis

Root Mean Square Fluctuation (RMSF) is a widely used metric to assess the flexibility of a system over time. It measures the average deviation of atomic positions from their equilibrium positions, providing insight into the mobility and stability of individual atoms or regions within a molecule. To compare the impact of Pan and 13a on HDAC3, we assessed the RMSF of Cα atoms.

As shown in Fig. 3A, the RMSF of the apo HDAC3 system differs from those of the HDAC3–Pan and HDAC3–13a systems. Overall, most residues in apo HDAC3 exhibit lower dynamicity, indicating a relatively stable structure, with RMSF values ranging from a minimum of 1.69 Å (residue 132) to a maximum of 15.28 Å (residue 398). However, elevated RMSF values are observed in the 71–101 and 391–426 regions, reflecting greater structural flexibility. These findings suggest that these regions may serve as flexible segments potentially involved in conformational changes, ligand binding, or interactions with regulatory proteins.

Fig. 3 (A) The RMSF diagrams of HADC3, HDAC3-17, HDAC3-11h. (B) The secondary structure changes of HADC3, HDAC3–Pan, HDAC3–13a in residues 71–101. (C) The secondary structure changes of HADC3, HDAC3–Pan, HDAC3–13a in residues 391–426. Docking results of Pan. The amino acid sequence of the nuclear receptor corepressor 2 (NCoR2) chain was renumbered to follow HDAC3 (residues 2–370), with residues 408–476 becoming residues 371–439.

In the HDAC3–Pan system, the RMSF values are generally higher than those in the apo HDAC3 system, particularly in the 351–401 region, where the maximum value reaches 24.23 Å at residue 398. This indicates that the binding of Pan increases the dynamicity and conformational changes of HDAC3, particularly in the 351–401 region. Additionally, the 71–101 region of the HDAC3–Pan system also shows higher fluctuations compared to the apo HDAC3, which may be related to the binding of Pan and may further enhance structural flexibility in this region.

The HDAC3–13a system also exhibits noticeable fluctuations in several key regions, especially the 351–401 region, where the maximum value reaches 28.12 Å at residue 398, which is significantly higher than the corresponding value in the apo HDAC3 system. Similar to Pan, 13a also leads to higher fluctuations in the 71–101 region. This suggests that the binding of the 13a ligand induces greater conformational changes and increased dynamicity in these regions.

Define Secondary Structure of Proteins (DSSP) analysis is a method for determining the secondary structure of proteins by evaluating the geometric properties of backbone hydrogen bonds. Residues are classified into different secondary structure types based on these calculations. In Fig. 3B and C, different colors are used to represent various structural elements. In the DSSP analysis of the apo HDAC3, HDAC3–Pan, and HDAC3–13a systems, we focused on secondary structure changes in the 71–101 region and the 391–426 region.

As shown in Fig. 3B, in the apo HDAC3 system, the 79–87 region primarily adopts a relatively stable α-helix structure. After the binding of the Pan, disruption of the helical structure is observed, particularly around residues 78–80, where a 3–10 helix primarily forms, connecting to the α-helix spanning residues 82–87. Similarly, the binding of the 13a ligand induces comparable structural changes, with the formation of a 3–10 helix linking to the adjacent α-helix. Additionally, in the 91–94 region, the apo HDAC3 system adopts a relatively stable anti-parallel conformation between 350 and 400 ns. However, in both the Pan and 13a systems, no long-lasting anti-parallel conformation is observed.

As shown in Fig. 3C, in the 391–426 region, all three systems exhibit similar secondary structures, with only a few areas showing differences. In the apo HDAC3 system, the 422–423 region transitions from a relatively stable bend conformation to a more stable turn structure after 230 ns. In both the HDAC3–Pan and HDAC3–13a systems, the 422–423 region predominantly consists of a turn structure. Pan and 13a induce similar structural alterations in HDAC3. Additionally, in the HDAC3–13a system, unlike in the other systems, the α-helix near residue 391 is disrupted, leading to a long-lasting bend structure.

Dynamic cross-correlation matrix analysis

Dynamic Cross-Correlation Matrix (DCCM) analysis is a powerful method for investigating the internal dynamic correlations of protein residues. By computing the displacement correlations between different residues, DCCM provides insights into the cooperative or anticorrelated motions within a protein. The correlation values range from −1 to 1, where values approaching 1 indicate a strong positive correlation (coordinated motion), values near −1 suggest strong anticorrelation (motions in opposite directions), and values close to 0 imply little to no correlation, signifying independent movements.

Through DCCM analysis of the apo HDAC3, HDAC3–Pan, and HDAC3–13a systems (Fig. 4), we can assess how ligand binding influences the protein's internal dynamics.

Fig. 4 The dynamical cross-correlation matrix diagrams of (A) HADC3 and (B) HDAC3–Pan, (C) HDAC3–13a. The amino acid sequence of the NCoR2 chain was renumbered to follow HDAC3 (residues 2–370), with residues 408–476 becoming residues 371–439.

In the apo HDAC3 system, a certain degree of positive correlation is observed along the diagonal, while correlations in regions further from the diagonal are weaker. Additionally, several regions of negative correlation are observed between the NCoR2 chain (residues ≥ 371) and the HDAC3 (residues < 371) chain. This observed negative correlation suggests the presence of inter-chain dynamic coupling. Such anti-correlated motions may reflect allosteric communication between the two proteins, where conformational changes in NCoR2 could influence the structural dynamics or functional state of HDAC3, possibly contributing to the regulation of its enzymatic activity.

In the HDAC3–Pan system, the positive correlations in regions further from the diagonal are weaker compared to the apo HDAC3 system, suggesting that ligand binding dampens long-range coordinated motions. This reduction in distal correlated dynamics implies that the presence of Pan stabilizes the overall structure by restricting large-scale conformational flexibility. In addition, the reduced extent of negative correlation between the NCoR2 and HDAC3 chains in the HDAC3–Pan system implies a weakened inter-chain dynamic coupling, suggesting that ligand binding reduces the allosteric communication between the two proteins.

Similar to the HDAC3–Pan system, the HDAC3–13a system also exhibits a weaker negative correlation between the NCoR2 and HDAC3 chains compared to the apo HDAC3 system. This suggests that 13a, like Pan, reduces the allosteric communication between the two proteins. However, unlike Pan, the binding of 13a does not lead to weaker positive correlations in regions further from the diagonal. Instead, the HDAC3–13a system retains positive correlations between residues within the HDAC3 chain (residues < 371) that are comparable to those observed in the apo HDAC3 system, with slightly stronger correlations in regions 100–150.

Distance analysis

The carbonyl oxygen atoms of the inhibitors Pan and 13a may coordinate with the zinc ion in the catalytic pocket, potentially influencing their binding interactions with the protein. To assess the likelihood of such coordination, we monitored the distances between the carbonyl oxygen atoms of both inhibitors and the catalytic zinc ion throughout the 500 ns GaMD simulation.

As shown in Fig. 5, the distance between the carbonyl oxygen atom of Pan and the zinc ion in the HDAC3 catalytic pocket exhibited greater fluctuations throughout the simulation. During the initial phase, the distance typically ranged from 2 to 5 Å, which is relatively long and suggests an unstable or weak coordination interaction. As the simulation progressed, the distance became somewhat more stabilized, mostly fluctuating between 2 and 4 Å. This range implies that Pan may form a transient or weak coordination bond with the zinc ion. However, the persistent fluctuations suggest that the interaction is not stably maintained throughout the simulation, possibly reflecting dynamic binding behavior or suboptimal positioning of the coordinating group.

Fig. 5 (A) Distance between carbonyl oxygen of ligand and zinc ion during the 500 ns simulations for HDAC3–Pan and HDAC3–13a. (B) The representative structure of HDAC3–Pan. (C) The representative structure of HDAC3–13a.

In contrast, the distance between the carbonyl oxygen atom of the 13a inhibitor and the zinc ion in the HDAC3 catalytic pocket remained relatively stable throughout the simulation. During the first 150 ns, the distance mostly fluctuated between 2 and 4 Å. As the simulation progressed, it gradually converged to a narrower range around 2–2.5 Å. This consistent proximity suggests a favorable positioning of the carbonyl group for stable coordination with the zinc ion. The stable and short distance implies that 13a is more likely to form and maintain a strong coordination bond with the zinc ion, potentially contributing to its enhanced inhibitory effect on HDAC3 compared to Pan.

PCA analysis

To determine how ligand binding affects the conformational distribution of HDAC3 and to identify the lowest-energy conformations in the HDAC3–Pan and HDAC3–13a systems, principal component analysis (PCA) and free energy landscape (FEL) construction were performed.

As shown in Fig. 6, in the HDAC3–Pan system, the low-energy states are structurally similar and can be categorized as a single cluster. We extracted the representative structure (lowest-energy conformation) from this cluster for further analysis. Consistent with our previous secondary structure analysis, the originally continuous α-helix in the 78–83 region is disrupted, forming two connected helical segments. Additionally, residues 175–180 exhibit a transition between an α-helix and a random coil, while residues 208–213 adopt a well-defined β-sheet conformation.

Fig. 6 Free energy landscape (FEL) and representative conformation of HDAC3–Pan system. The amino acid sequence of the NCoR2 chain was renumbered to follow HDAC3 (residues 2–370), with residues 408–476 becoming residues 371–439.

Fig. S6† shows the evolution of interactions between Pan and the surrounding residues of HDAC3 over the course of the simulation. Fig. S7† illustrates the interactions observed in the representative structure.

Pan forms hydrogen bonds with HIS135, GLY143, and TYR298, pi–pi interactions with PHE144 and PHE200, and a pi–alkyl interaction with PRO23. Additionally, it engages in van der Waals interactions with nearby residues such as ASP92, ASP93, ASP170, and LEU266.

Compared to the docking results, the low-energy conformation of the MD simulation reveals both preserved and altered interactions. The pi–pi interactions with PHE144 and PHE200 and the pi–alkyl interaction with PRO23 remains intact, suggesting that these key interactions play a dominant role in Pan binding to HDAC3. However, GLY21 and ASP93, which initially formed a hydrogen bond in the docking model, instead engages in van der Waals interactions. Similarly, the hydrogen bond between Pan and HIS134 is lost, while new hydrogen bonds form with GLY143 and TYR298. This indicates that Pan undergoes subtle conformational adjustments within the binding pocket while maintaining hydrogen bond interactions with surrounding residues. Notably, in the low-energy conformation, no metal–acceptor interaction is detected between the zinc ion and the carbonyl oxygen atom, likely due to their relatively large separation distance, which is consistent with our previous distance analysis.

As shown in Fig. 7, in the HDAC3–13a system, the low-energy states can be categorized into three clusters. We extracted representative structures (the lowest-energy conformations) from each cluster for further analysis. Unlike the HDAC3–Pan system, all three representative conformations exhibit a well-preserved α-helix in the 175–180 region. Notably, in State A, the β-sheet spanning residues 218–233 is disrupted, whereas this structural alteration is not observed in State B or State C. Additionally, in State C, the 78–83 region adopts a combination of random coil and α-helix, similar to what is observed in the HDAC3–Pan system, whereas in States A and B, this segment maintains a continuous helical structure.

Fig. 7 Free energy landscape (FEL) and representative conformation of HDAC3–13a system. The lowest-energy conformation of each cluster is shown in (A)–(C), respectively. The amino acid sequence of the NCoR2 chain was renumbered to follow HDAC3 (residues 2–370), with residues 408–476 becoming residues 371–439.

Fig. S8† shows the evolution of interactions between 13a and the surrounding residues of HDAC3 over the course of the simulation. Fig. S9–S11† illustrates the interactions observed in the representative structures.

In State A (Fig. S9†), 13a forms hydrogen bonds with HIS135, ASP170, ASP259, and TYR298, pi–pi interactions with PHE144 and PHE200, and pi–alkyl interactions with PRO23, MET24, LEU133, and CYS145. Additionally, it engages in attractive charge interactions with ASP170 and ASP259, as well as a metal–acceptor interaction with the zinc ion. It also forms van der Waals interactions with nearby residues such as ASP92, ASP93, GLY132, and LEU266.

In State B (Fig. S10†), 13a forms hydrogen bonds with HIS135, ASP170, and TYR298, pi–pi interactions with PHE144 and PHE200, and pi–alkyl interactions with MET24, LEU133, PHE144, and CYS145. It also engages in an attractive charge interaction with ASP170 and a metal–acceptor interaction with the zinc ion. Additionally, van der Waals interactions are observed with ASP93, ASP259, GLY132, LEU266, etc.

In State C (Fig. S11†), 13a forms hydrogen bonds with ASP170 and TYR298, pi–pi interactions with PHE144 and PHE200, and pi–alkyl interactions with MET24, LEU133, CYS145, and LEU266. It also engages in an attractive charge interaction with ASP170 and a metal–acceptor interaction with the zinc ion. van der Waals interactions are observed with ASP93, ASP259, GLY143, GLY296, etc.

Compared to the docking results, the three representative low-energy conformations of the HDAC3–13a system exhibit both preserved and altered interactions. The pi–pi interactions between 13a and PHE144 remain intact across all three conformations. Additionally, the hydrogen bond with HIS135 is preserved in two of the three representative conformations, suggesting that these key interactions play a dominant role in 13a binding to HDAC3. However, other hydrogen bond interactions undergo rearrangements in the low-energy conformations. Notably, in all three representative structures, 13a forms hydrogen bonds with ASP170 and TYR298. Furthermore, 13a engages in attractive charge interactions with ASP170, highlighting its strong electrostatic contributions to binding. Additionally, the interaction between 13a and the zinc ion shifts from a van der Waals interaction in the docking model to a metal–acceptor interaction in the low-energy conformations, which aligns with our previous distance analysis. Throughout the HDAC3–13a system simulation, the distance between the zinc ion and the carbonyl oxygen atom remains relatively short and exhibits lower fluctuations compared to the HDAC3–Pan system.

Alanine mutation analysis

In this analysis, we mutated the residues surrounding the inhibitors (Pan and 13a) to alanine, effectively replacing active residues with alanine and removing the active functional groups from the side chains, substituting them with methyl groups that have a relatively minor impact on the protein structure. The objective was to investigate the influence of these active residues on the structural stability of the protein. The test results are presented in Tables 1 and 2.

Table 1 Alanine mutation of HDAC3–Pan

Mutation	Mutation energy (kcal mol^−1)	Effect
GLY21	−0.05	Neutral
HIS22	0.5	Neutral
PRO23	0.52	Destabilizing
MET24	−0.12	Neutral
ASP92	0.33	Neutral
ASP93	0.23	Neutral
CYS94	−0.01	Neutral
PRO95	0.04	Neutral
LEU133	−0.04	Neutral
HIS134	−0.88	Stabilizing
HIS135	0.2	Neutral
GLY143	−0.56	Stabilizing
PHE144	1.75	Destabilizing
PHE200	1.11	Destabilizing
LEU266	0.02	Neutral
GLY295	−0.05	Neutral
GLY296	0.77	Destabilizing
GLY297	−0.08	Neutral
TYR298	0.2	Neutral

---

Table 2 Alanine mutation of HDAC3–13a

Mutation	Mutation energy (kcal mol^−1)	Effect
GLY21	−0.05	Neutral
HIS22	0.14	Neutral
PRO23	1.06	Destabilizing
MET24	0.3	Neutral
ARG28	−0.68	Stabilizing
ASP92	0.27	Neutral
ASP93	0.23	Neutral
GLY132	0.08	Neutral
LEU133	1.02	Destabilizing
HIS134	−0.03	Neutral
HIS135	1.76	Destabilizing
GLY143	−0.03	Neutral
PHE144	1.76	Destabilizing
CYS145	0.03	Neutral
ILE150	−0.15	Neutral
PHE200	1.07	Destabilizing
GLN255	0.62	Destabilizing
ARG265	−0.26	Neutral
LEU266	1.17	Destabilizing
GLY295	1.05	Destabilizing
GLY296	0.1	Neutral
GLY297	0.05	Neutral
TYR298	−1.31	Stabilizing

---

In the HDAC3–Pan system, most of mutations showed minimal impact on the overall binding energy (mutation energy between −0.5 and 0.5 kcal mol⁻¹), suggesting that these residues are not essential for maintaining the stability of the complex. Notably, three mutations—PHE144 (1.75 kcal mol⁻¹), PHE200 (1.11 kcal mol⁻¹), and GLY296 (0.77 kcal mol⁻¹)—were classified as destabilizing, indicating that these residues may play key roles in maintaining local structural integrity or facilitating ligand interaction. The relatively high mutation energies suggest that substituting these residues with alanine disrupts favorable interactions, possibly through the loss of hydrophobic contacts or structural constraints. Conversely, mutations at HIS134 (−0.88 kcal mol⁻¹) and GLY143 (−0.56 kcal mol⁻¹) resulted in negative mutation energies, indicating a stabilizing effect. This may suggest that these residues contribute to local flexibility or unfavorable steric hindrance in the wild-type structure, which is alleviated upon alanine substitution.

In the HDAC3–13a system, mutations at PRO23, LEU133, HIS135, PHE144, PHE200, GLN255, LEU266, and GLY295 were predicted to be destabilizing, each showing mutation energies greater than 1.0 kcal mol⁻¹. These results suggest that these residues play critical roles in maintaining the structural integrity of the protein–ligand complex. In particular, HIS135 and PHE144 each exhibited a mutation energy of 1.76 kcal mol⁻¹, indicating that they are essential for stable binding with the 13a inhibitor. Conversely, ARG28 and TYR298 mutations were predicted to be stabilizing, with mutation energies of −0.68 and −1.31 kcal mol⁻¹, respectively. The stabilizing effect of the TYR298 mutation suggests it may help anchor the inhibitor in a favorable binding orientation, potentially enhancing interaction specificity or strength.

Discussion

This study provides a comprehensive comparison of the dynamic and structural consequences of HDAC3 inhibition by two compounds: Pan and 13a. Through a range of simulation analyses, we demonstrate that although both ligands target the same protein, they influence its behavior in notably different ways.

Structural stability: RMSD analyses

RMSD analysis reveals that Pan induces moderate conformational changes early in the simulation, after which the HDAC3–Pan complex reaches a relatively stable state. In contrast, 13a binding results in more sustained and pronounced fluctuations, suggesting a higher degree of flexibility or structural rearrangement. This observation is consistent with 13a's known slow-on/slow-off binding kinetics and implies a distinct mode of action compared to Pan.

Structural compactness: R_g and SASA analyses

Both inhibitors lead to a reduction in the R_g and SASA, indicating increased structural compactness of the HDAC3 complex. However, the HDAC3–13a complex undergoes a more marked compaction than the HDAC3–Pan complex. This suggests that 13a binding triggers certain regions of the protein to close or undergo structural rearrangement. This compactification of the structure may contribute to the protein's stability and functionality in specific biological environments.

Flexibility and local conformational dynamics: RMSF and DSSP analyses

RMSF analysis shows increased flexibility in both ligand–bound complexes compared to apo HDAC3. Among them, the HDAC3–13a complex exhibits the highest degree of fluctuation, indicating a more dynamic interaction between HDAC3 and 13a. This heightened flexibility may be associated with its slow-on/slow-off binding kinetics and distinct inhibitory mechanism.

Furthermore, only the 13a–bound complex exhibits a persistent disruption of the α-helix near residue 391, leading to a stable bent structure. This conformational change may be linked to the therapeutic mechanism of 13a, as the loss of α-helical integrity could affect HDAC3's interactions with other proteins, potentially influencing its therapeutic efficacy in AML treatment.

Allosteric and intra-protein communication: DCCM analysis

Dynamical cross-correlation matrix (DCCM) analysis shows that both inhibitors disrupt negative inter-chain correlations between HDAC3 and NCoR2, indicating interference with allosteric communication. Interestingly, Pan binding also dampens long-range intra-protein correlations within HDAC3, suggesting a rigidification of the protein structure that may lock it in an inactive conformation. In contrast, 13a maintains the large-scale flexibility characteristic of apo HDAC3, reflecting its potentially unique regulatory mechanism.

Zinc ion coordination: distance analysis

Distance measurements between the inhibitors and the catalytic zinc ion highlight that 13a maintains a consistently short and stable distance, indicative of a robust coordination bond. Although Pan sometimes came into close proximity with the zinc ion, its interaction remained more variable and less stable overall. This difference in binding stability likely contributes to 13a's greater inhibitory potency.

Principal component analysis (PCA) and binding mode

Unlike Pan, which stabilizes HDAC3 into relatively uniform low-energy conformations, 13a binding leads to multiple distinct conformational states. This observation aligns with the RMSD results and reflects the structural complexity introduced by 13a. The presence of several low-energy conformations and the potential transitions among them may explain why 13a behaves as a slow-on/slow-off inhibitor.

Analysis of low-energy structures further shows that both Pan and 13a induce similar secondary structural changes in HDAC3 to some extent, but key differences are also evident. In particular, 13a appears to promote more substantial conformational rearrangements upon binding—such as those observed in Fig. 7A—which may contribute to its unique inhibitory mechanism.

In addition, in the low-energy conformations of both HDAC3–13a and HDAC3–Pan complexes, the inhibitors form key interactions with residues such as HIS135, PHE144, PHE200, and TYR298, highlighting the importance of these residues for ligand binding. Compared to Pan, 13a forms more π–alkyl interactions, along with additional attractive charge and metal–acceptor interactions. These enhanced interactions may underlie 13a's superior inhibitory potency against HDAC3.

Binding energy contributions: alanine scanning mutagenesis

Alanine scanning underscores the differential binding dependencies of the two inhibitors. The HDAC3–Pan complex appears less affected by single-point mutations, indicating a more tolerant and flexible interaction network. In contrast, the HDAC3–13a complex is highly sensitive to mutations at specific residues, notably HIS135, PHE144, and PHE200, reinforcing the idea that 13a relies on more precise and critical contacts for stable binding.

Concluding Remarks

Collectively, these results highlight how distinct inhibitors modulate HDAC3 dynamics and structure through different mechanisms. Pan tends to stabilize HDAC3 in a more rigid, less flexible conformation, potentially favoring immediate but transient inhibition. Conversely, 13a induces significant conformational changes, enhances protein compactness, and maintains dynamic flexibility, characteristics that may underlie its slow-on/slow-off kinetics and stronger inhibitory performance. These insights offer valuable implications for the rational design of HDAC3-targeted therapeutics, particularly in the context of diseases such as AML where selective and sustained inhibition is desired.

Materials and methods

System preparation

The 3D structure of HDAC3 (PDB: 4A69) was obtained from the Protein Data Bank.³⁰ Since 4A69 contains two identical HDAC3 chains, one of them was removed during protein preparation in UCSF Chimera.³⁴ The NCoR2 chain, which interacts with the retained HDAC3 chain, was preserved. Additionally, all ligands and water molecules were removed. The structures of Panobinostat and compound 13a were modeled and optimized using Discovery Studio.³⁵ AutoDock Vina v1.2.6 (ref. 31) was then employed to dock Panobinostat and 13a to HDAC3, resulting in the HDAC3–Pan and HDAC3–13a complexes. The docking grid was defined based on the active site information from the HDAC3 crystal structure (PDB ID: 4A69), including the zinc ion and surrounding key residues: HIS134, HIS135, GLY143, ASP170, HIS172, ASP259, GLY296, and TYR298. Vina simplifies the process by eliminating the need to manually select atomic types and predict grid maps for ligands and receptors, thereby improving the overall accuracy of binding predictions and accelerating the search process with more straightforward scoring functions. The structure with the lowest binding affinity from the docking results was selected as the starting point for molecular dynamics (MD) simulations.

To parameterize the Zn²⁺ coordination center in HDAC3, we used the MCPB.py³⁶ module provided by AMBER. The metal site and its directly coordinated residues/ligands were extracted from the crystal structure to build a small model. This model was then optimized using Gaussian 16 at the B3LYP/6-31G* level of theory. Following geometry optimization, electrostatic potential (ESP) data were generated and RESP charges were fitted. Using these results, MCPB.py produced the corresponding force field parameters (frcmod) and library (mol2/lib) files. The final metal center model was integrated into the full protein topology for subsequent molecular dynamics simulations. The ligand was treated in the same manner, with geometry optimization and RESP charge fitting performed using Gaussian.

The protein–inhibitor complex was prepared using LEaP module in AMBER 22.³⁷ The ff19SB³⁸ force field was applied to the protein, and the ligand was parameterized using the GAFF2 force field. The Zn²⁺ metal center and its coordinating residues were treated using parameters generated from the MCPB.py module. The resulting complex was solvated in a truncated octahedral box of OPC water³⁹ molecules with a 12 Å buffer. To neutralize the net charge of the system, 6 Na⁺ counterions were added. The final system was composed of the protein, the metal center, the bound ligand, solvent molecules, and counterions. Three systems were constructed: HDAC3 (86 206 atoms), HDAC3–Pan (80 875 atoms), and HDAC3–13a (80 853 atoms).

Equilibrium simulation

Prior to the GaMD production runs, a series of conventional molecular dynamics simulations were performed using the AMBER 22 simulation package. This process was taken to prevent incorrect positioning of certain atoms from docking and to provide a more reasonable and stable starting structure for the subsequent GaMD simulations. The protocol consisted of the following steps:

Minimization. The system was subjected to energy minimization consisting of a total of 10 000 cycles, including 5000 steps of the steepest descent algorithm followed by 5000 steps of the conjugate gradient algorithm. This step aimed to remove unfavorable contacts and relax the initial structure.

Heating. Velocities were assigned to all atoms according to a Boltzmann distribution, and the system was gradually heated from 0.1 K to 310 K using a linear temperature ramp. The heating phase was conducted under the NVT ensemble for 25 000 steps with a 2 fs time step, corresponding to a total simulation time of 50 ps. Langevin dynamics with a collision frequency of 2.0 ps⁻¹ was used to regulate temperature. The gradual temperature increase was designed to ensure a smooth and stable transition to the target temperature.

Density equilibration. Following the heating step, the system was equilibrated under isothermal-isobaric (NPT) conditions for 50 ps to stabilize the density. This phase consisted of 25 000 steps with a 2 fs time step at a constant temperature of 310 K using Langevin dynamics (collision frequency = 2.0 ps⁻¹). Pressure was maintained at 1 atm using an isotropic Berendsen barostat with a relaxation time of 1.0 ps.

Throughout the minimization, heating, and density equilibration steps, positional restraints were applied to all protein atoms using harmonic potentials with a force constant of 2.0 kcal mol⁻¹·Å⁻² to preserve the structural integrity of the protein.

Final equilibration. After the initial density equilibration, a final equilibration phase was conducted under NPT conditions for 500 ps (250 000 steps, 2 fs time step) to ensure the stability of temperature, pressure, and volume. Langevin dynamics with a collision frequency of 2.0 ps⁻¹ was employed for temperature control, and pressure was maintained at 1 atm using an isotropic Berendsen barostat with a relaxation time of 2.0 ps. During this stage, positional restraints were removed, allowing the protein and solvent to equilibrate freely within a thermodynamically relevant ensemble.

Production MD simulation. Following equilibration, the production MD simulations were performed under NPT conditions at 310 K and 1 atm pressure. The simulations were conducted using a 2 fs time step. Langevin dynamics with a collision frequency of 2.0 ps⁻¹ was employed for temperature regulation, while pressure was controlled using a Monte Carlo barostat with a relaxation time of 5.0 ps. Periodic boundary conditions were applied to mimic bulk solvent behavior.

All simulations were conducted using the pmemd.cuda module in AMBER for GPU-accelerated performance. A nonbonded cutoff of 8.0 Å was applied, and long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME⁴⁰) method. The SHAKE⁴¹ algorithm was used to constrain all covalent bonds involving hydrogen atoms, allowing the use of a 2 fs integration time step.

Gaussian accelerated molecular dynamics simulation

Gaussian accelerated molecular dynamics (GaMD)⁴² enhances the conformational sampling of biomolecules. It works by adding a harmonic boost potential to smooth the biomolecular potential energy surface and reduce energy barriers, thereby accelerating transitions between different low-energy states. For a system with N atoms located at r⃑ = {r⃑₁,⋯r⃑_N},when the system potential V(r⃑) is below the threshold energy E, a lifting potential is added as:

Among them, k is the harmonic force constant. The revised system potential is:

when the system potential energy is higher than the threshold value, that is, Vr ≥ E, the boost potential is set to zero, in which case V*(r⃑) = V(r⃑).

In order to smooth the potential energy surface of enhanced sampling, the following conditions need to be met: First, for any two potential values V₁(r⃑) and V₂(r⃑) on the original energy surface, if V₁(r⃑) < V₂(r⃑), ΔV should be a monotone function that does not change the relative order of the bias potential values. By bringing V*(r⃑) into (2) we get:

Second, if V₁(r⃑) < V₂(r⃑), the potential difference observed on a smooth energy surface should be smaller than the original potential difference. Similarly, bringing V*(r⃑) into (2), we can derive:

We need to combine eqn (3) and (4) to set the threshold energy E to the following range:

Among them, V_min and V_max are the minimum and maximum potential energy of the system, respectively. In order for eqn (5) to be valid, k needs to satisfy:

k₀ is defined to determine the magnitude of the applied boost potential. Specifically, k ≡ k₀(1/(V_max − V_min)); then 0 < k₀ ≤ 1.

Third, the standard deviation of ΔV needs to be small enough (i.e. narrow distribution) to ensure accurate reweighting using second-order cumulative extension:

Among them, V_avg and σ_V are the mean and standard deviation of the potential energy of the system, σ_ΔV is the standard deviation of ΔV, and σ₀ is the upper limit specified by the user so that the weights can be adjusted accurately.

The setting formula (5) gives the value range of the threshold energy E. When E is set to E = V_max, E and k are substituted to get:

The larger σ_ΔV is obtained from the original potential energy surface (especially for large biomolecules). When the threshold energy E = V_min + 1/k is set as its upper bound according to eqn (5) n, E and k are substituted into eqn (7), thus:

From given E and k₀, we can calculate the boost potential as:

Similar to aMD, GaMD offers the option to add only the total potential boost ΔV_p, only the dihedral potential boost ΔV_D, or the bipotential boost (ΔV_p and ΔV_D). For enhanced sampling, dual-boost simulations generally provide higher acceleration than the other two types of simulations. The simulation parameters are composed of threshold energy value and effective harmonic force constants k_0P and k_0D.

Analysis of MD trajectories

All trajectory analyses were performed using the cpptraj module of the AMBER software suite.

Principal component analysis (PCA) was performed to capture the major conformational motions sampled during accelerated molecular dynamics (aMD) simulations. The PCA was performed on the atomic Cartesian coordinates of residues 1–438 (covering the full protein) excluding hydrogens, as extracted from the GaMD trajectory. Prior to PCA, the trajectory was aligned to the first frame, and the average structure was used for covariance matrix calculation. The covariance matrix was subsequently diagonalized to obtain eigenvectors (principal components, PCs) and their corresponding eigenvalues, which describe the dominant modes of collective motion. The first two principal components (PC1 and PC2) were selected as reaction coordinates for reconstructing the conformational free energy surfaces, as they capture the dominant collective motions and account for a significant proportion of the total variance—30.33% for the HDAC3–Pan system and 40.08% for the HDAC3–13a system.

To correct for the bias introduced by aMD and recover the canonical ensemble distribution, energetic reweighting was applied using the PyReweighting toolkit⁴³ (Miao et al.). In this study, the second-order cumulant expansion (CE2) method was used. The boost potential values extracted from the aMD simulation log were combined with the PC1–PC2 projections to calculate a reweighted two-dimensional free energy surface. The resulting landscape provides insights into the conformational preferences and energetics of the system under study.

We used Discovery Studio and the Python packages MDAnalysis⁴⁴ and ProLIF⁴⁵ to analyze protein–small molecule interactions during the MD simulations.

Alanine scanning

Using Discovery Studio 2021, we conducted an alanine scanning mutagenesis study to investigate the active sites of both the HDAC3–Pan and HDAC3–13a complexes. This technique was employed to systematically replace specific amino acids near the active sites with alanine. The rationale behind this approach is that alanine, being a small and non-reactive amino acid, can be used to evaluate the importance of side-chain interactions in inhibitor binding. By replacing key residues with alanine, we aimed to identify which amino acids play a critical role in stabilizing the inhibitor within the binding pocket. This process allowed us to pinpoint the residues that contribute significantly to the binding affinity and overall inhibitory effect of the compounds on HDAC3. The results from the alanine scanning provided valuable insights into the molecular interactions that are crucial for inhibitor binding, thereby highlighting potential targets for enhancing the efficacy of HDAC3 inhibitors.

The mutation energy (ΔΔG_mut), which reflects the impact of each mutation on binding affinity, is determined by subtracting the binding free energy of the wild-type protein from that of the mutated structure:

ΔΔG_mut = ΔG_bind(mutant) − ΔG_bind(wild type) (12)

The binding free energy, ΔG_bind, represents the difference between the free energy of the complex and the unbound state. All energy components are evaluated using CHARMm, with the electrostatic energy computed via a Generalized Born implicit solvent model. The overall energy is derived as an empirically weighted combination of van der Waals interactions (E_vdW), electrostatic contributions (ΔG_elec), an entropy term (–TS_sc) associated with alterations in side-chain flexibility, and a non-polar, surface-area-dependent solvation energy component (ΔG_np).

Conclusions

This study utilized Gaussian accelerated molecular dynamics simulations to compare the binding mechanisms of Panobinostat and compound 13a with HDAC3, a key target in cancer therapy. Our results demonstrate that 13a exhibits superior inhibitory activity due to its stronger and more stable interactions with HDAC3, including enhanced coordination with the catalytic zinc ion and a greater number of hydrogen bonds and electrostatic interactions. Notably, 13a induces significant conformational flexibility in HDAC3, stabilizing the enzyme in multiple distinct states, which likely contributes to its slow-on/slow-off binding kinetics. These findings provide a molecular basis for the enhanced potency of 13a and underscore its potential as a therapeutic agent for acute myelogenous leukemia (AML). Moreover, the identification of critical residues through alanine scanning offers valuable insights for the design of more effective HDAC3 inhibitors. This work highlights the utility of advanced computational methods in elucidating enzyme-inhibitor dynamics and paves the way for the development of next-generation cancer therapeutics.

Data availability

All data generated or analysed during this study are included in this published article.

Author contributions

Xia Yu: contextualization, experimentation & manuscript drafting. Hengzheng Yang: experimentation, data collection, interpretation, analysis & manuscript drafting. Baiji Xue: manuscript drafting, data analysis & reviewing. Tong Liu: manuscript drafting & editing. Xue Zhang: manuscript drafting & experimental supervision. Yang Xu: manuscript drafting, data analysis & reviewing. Xueliang Zhao: manuscript drafting & editing. Xianwen Yue: manuscript drafting, editing & supervision.

Conflicts of interest

The authors declare no competing interests.

Acknowledgements

Thank you to DA Case for providing Amber22 for free. This research was funded by Jilin Provincial Science and Technology Department (No. YDZJ202301ZYTS19) & The key research and development project of Baicheng Glycyrrhiza of Baicheng Medical College (No. BCGC202413).

References

1. A. G. Muntean and J. L. Hess, Am. J. Pathol., 2009, 175, 1353–1361.
2. Y. Cheng, C. He, M. Wang, X. Ma, F. Mo, S. Yang, J. Han and X. Wei, Signal Transduction Targeted Ther., 2019, 4, 62.
3. S. Ropero and M. Esteller, Mol. Oncol., 2007, 1, 19–25.
4. M. Haberland, R. L. Montgomery and E. N. Olson, Nat. Rev. Genet., 2009, 10, 32–42.
5. J. Juárez-Jiménez, A. A. Gupta, G. Karunanithy, A. S. J. S. Mey, C. Georgiou, H. Ioannidis, A. De Simone, P. N. Barlow, A. N. Hulme, M. D. Walkinshaw, A. J. Baldwin and J. Michel, Chem. Sci., 2020, 11, 2670–2680.
6. D.-Q. Chen, C. Yu, X.-F. Zhang, Z.-F. Liu, R. Wang, M. Jiang, H. Chen, F. Yan, M. Tao, L.-B. Chen, H. Zhu and J.-F. Feng, Ther. Adv. Med. Oncol., 2018, 10, 1758835918783132.
7. N. Wu, Q. Sun, L. Yang, H. Sun, Z. Zhou, Q. Hu, C. Li, D. Wang, L. Zhang, Y. Hu and X. Cong, Int. J. Biol. Macromol., 2025, 300, 140310.
8. H.-Y. Hsieh, H.-C. Chuang, F.-H. Shen, K. Detroja, L.-W. Hsin and C.-S. Chen, Eur. J. Med. Chem., 2017, 140, 42–51.
9. H. Wang, Y.-C. Liu, C.-Y. Zhu, F. Yan, M.-Z. Wang, X.-S. Chen, X.-K. Wang, B.-X. Pang, Y.-H. Li, D.-H. Liu, C.-J. Gao, S.-J. Liu and L.-P. Dou, J. Exp. Clin. Cancer Res., 2020, 39, 278.
10. Y. Li, Y. Liu, N. Zhao, X. Yang, Y. Li, F. Zhai, X. Zang and W. Cui, Cell Death Dis., 2020, 11, 753.
11. R. Sarkar, S. Banerjee, S. A. Amin, N. Adhikari and T. Jha, Eur. J. Med. Chem., 2020, 192, 112171.
12. L. Zhang, Y. Chen, Q. Jiang, W. Song and L. Zhang, Eur. J. Med. Chem., 2019, 162, 534–542.
13. R. L. Siegel, K. D. Miller and A. Jemal, Ca-Cancer J. Clin., 2020, 70, 7–30.
14. P. Mehdipour, F. Santoro, O. A. Botrugno, M. Romanenghi, C. Pagliuca, G. M. Matthews, R. W. Johnstone and S. Minucci, Leukemia, 2017, 31, 995–997.
15. J. Long, W. Y. Fang, L. Chang, W. H. Gao, Y. Shen, M. Y. Jia, Y. X. Zhang, Y. Wang, H. B. Dou, W. J. Zhang, J. Zhu, A. B. Liang, J. M. Li and J. Hu, Leukemia, 2017, 31, 2761–2770.
16. X. Li, Y. Jiang, Y. K. Peterson, T. Xu, R. A. Himes, X. Luo, G. Yin, E. S. Inks, N. Dolloff, S. Halene, S. S. L. Chan and C. J. Chou, J. Med. Chem., 2020, 63, 5501–5525.
17. M. Duvic and J. Vu, Biologics, 2007, 1, 377–392.
18. E. M. Bertino and G. A. Otterson, Expert Opin. Invest. Drugs, 2011, 20, 1151–1158.
19. K. Ding, H. Liu, H. Yang, H. Zhu, J. Ma, H. Peng, H. Huang, W. Shi, L. Cao, W. Wu, X. Zhao, X. Shi, J. Li, X. Zhang and L. Fan, Med, 2024, 5, 1393–1401.e2.
20. S. L. Greig, Target. Oncol., 2016, 11, 107–114.
21. M. D. M. Traoré, V. Zwick, C. A. Simões-Pires, A. Nurisso, M. Issa, M. Cuendet, M. Maynadier, S. Wein, H. Vial, H. Jamet and Y.-S. Wong, ACS Omega, 2017, 2, 1550–1562.
22. C. Micelli and G. Rastelli, Drug Discovery Today, 2015, 20, 718–735.
23. Y. Wang, R. L. Stowe, C. E. Pinello, G. Tian, F. Madoux, D. Li, L. Y. Zhao, J.-L. Li, Y. Wang, Y. Wang, H. Ma, P. Hodder, W. R. Roush and D. Liao, Chem. Biol., 2015, 22, 273–284.
24. S. I. Son, J. Cao, C.-L. Zhu, S. P. Miller and H. Lin, ACS Chem. Biol., 2019, 14, 1393–1397.
25. X. Li, Y. Jiang, Y. K. Peterson, T. Xu, R. A. Himes, X. Luo, G. Yin, E. S. Inks, N. Dolloff, S. Halene, S. S. L. Chan and C. J. Chou, J. Med. Chem., 2020, 63, 5501–5525.
26. F. Guo, H. Yang, X. Bai, J. Li, W. Han and W. Li, J. Biomol. Struct. Dyn., 2024, 42, 13779–13792.
27. J. Chen, Q. Zeng, W. Wang, Q. Hu and H. Bao, RSC Adv., 2022, 12, 1742–1757.
28. H. N. Do, S. Akhter and Y. Miao, Front. Mol. Biosci., 2021, 8, 673170.
29. W. Yang, J. Wang, L. Zhao and J. Chen, Molecules, 2024, 29, 3377.
30. P. J. Watson, L. Fairall, G. M. Santos and J. W. R. Schwabe, Nature, 2012, 481, 335–340.
31. O. Trott and A. J. Olson, J. Comput. Chem., 2010, 31, 455–461.
32. H. Kumar and M. E. Sobhia, ACS Med. Chem. Lett., 2024, 15, 1306–1318.
33. D. K. Patel, H. Kumar and M. E. Sobhia, Comput. Biol. Chem., 2024, 112, 108132.
34. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng and T. E. Ferrin, J. Comput. Chem., 2004, 25, 1605–1612.
35. BIOVIA, Dassault Systèmes, Discovery Studio, Version 2021, San Diego: Dassault Systèmes, 2021.
36. P. Li and K. M. Merz, J. Chem. Inf. Model., 2016, 56, 599–604.
37. D. A. Case, H. M. Aktulga, K. Belfon, D. S. Cerutti, G. A. Cisneros, V. W. D. Cruzeiro, N. Forouzesh, T. J. Giese, A. W. Götz, H. Gohlke, S. Izadi, K. Kasavajhala, M. C. Kaymak, E. King, T. Kurtzman, T.-S. Lee, P. Li, J. Liu, T. Luchko, R. Luo, M. Manathunga, M. R. Machado, H. M. Nguyen, K. A. O'Hearn, A. V. Onufriev, F. Pan, S. Pantano, R. Qi, A. Rahnamoun, A. Risheh, S. Schott-Verdugo, A. Shajan, J. Swails, J. Wang, H. Wei, X. Wu, Y. Wu, S. Zhang, S. Zhao, Q. Zhu, T. E. I. Cheatham, D. R. Roe, A. Roitberg, C. Simmerling, D. M. York, M. C. Nagan and K. M. Jr. Merz, J. Chem. Inf. Model., 2023, 63, 6183–6191.
38. C. Tian, K. Kasavajhala, K. A. A. Belfon, L. Raguette, H. Huang, A. N. Migues, J. Bickel, Y. Wang, J. Pincay, Q. Wu and C. Simmerling, J. Chem. Theory Comput., 2020, 16, 528–552.
39. S. Izadi, R. Anandakrishnan and A. V. Onufriev, J. Phys. Chem. Lett., 2014, 5, 3863–3871.
40. M. J. Harvey and G. De Fabritiis, J. Chem. Theory Comput., 2009, 5, 2371–2377.
41. V. Kräutler, W. F. van Gunsteren and P. H. Hünenberger, J. Comput. Chem., 2001, 22, 501–508.
42. Y. Miao, V. A. Feher and J. A. McCammon, J. Chem. Theory Comput., 2015, 11, 3584–3595.
43. Y. Miao, W. Sinko, L. Pierce, D. Bucher, R. C. Walker and J. A. McCammon, J. Chem. Theory Comput., 2014, 10, 2677–2689.
44. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
45. C. Bouysset and S. Fiorucci, J. Cheminf., 2021, 13, 72.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ra01129a

‡ These authors contributed equally to this work.

---