Exploration of Natural Products for Targeting IDH1/2 Mutations in Acute Myeloid Leukemia Through Ligand-Based Pharmacophore Screening,Docking,ADME-T,and Molecular Dynamic Simulation Approaches

Abstract

Background:

Mutations in isocitrate dehydrogenase 1 (IDH1) and 2 (IDH2) are prevalent drivers of acute myeloid leukemia (AML). While targeted therapies exist, resistance can emerge. This study explored the potential of natural products to identify novel dual IDH inhibitors.

Methods:

In silico screening of the COCONUT database was performed using Lipinski’s Rule of Five. Pharmacophore modeling identified crucial features for IDH binding. Docking simulations with Glide (Schrödinger) assessed binding affinity, followed by MM-GBSA calculations for free energy estimation. The most promising candidate underwent ADME/T and toxicity analysis. Finally, molecular dynamics (MD) simulations evaluated the stability of protein-ligand complexes and binding interactions, followed by trajectory analysis using dynamical cross-correlation matrix (DCCM) and principal component analysis (PCA).

Results:

Ternstroside D (CNP0166496) emerged as a potential dual inhibitor of IDH1 and IDH2 mutations. Docking and MM-GBSA analyses showed strong affinities with IDH1 (–14.2, –84.45 kcal/mol) and IDH2 (–16.8, –60.73 kcal/mol), exceeding those of reference inhibitors GSK321A (–9.6 kcal/mol) and Enasidenib (–8.9 kcal/mol). Key hydrogen-bond interactions with catalytic residues and stable binding during MD simulations support its dual mechanism. ADME/T predictions indicated drug-like properties and a favorable safety profile, highlighting Ternstroside D as a natural scaffold with superior binding compared with existing IDH inhibitors.

Conclusion:

This in silico study provides compelling evidence for Ternstroside D (CNP0166496) as a promising dual inhibitor for IDH1 and IDH2 mutations in AML. Furthermore, in vitro and in vivo studies are warranted to validate these findings.

Introduction

Acute myeloid leukemia is an aggressive cancer affecting hematopoietic progenitor or stem cells, marked by an abnormal surge in myeloblast production. Consequently, this diminishes fully matured cellular components like erythrocytes, thrombocytes, and leukocytes, giving rise to critical health complications. ¹ In 2022, the United States was estimated to witness around 60 650 fresh instances of leukemia, as per predictions by the American Cancer Society.^2,3 Among its various forms, AML emerged as the deadliest, contributing to nearly 48% of fatalities linked to this blood cancer. ⁴ Although chemotherapy serves as the primary therapeutic approach, the prognosis for individuals with AML remains poor, with a survival probability of merely 28% over 5 years. In elderly patients or those with additional health conditions, the median lifespan ranges between 5 and 10 months. ⁵

Acute myeloid leukemia is often associated with the dysregulated expression of certain genes, including the isocitrate dehydrogenase (IDH) genes, which play a crucial role in the disease’s pathogenesis. Mutations in IDH1 and IDH2 were first identified in AML in 2008 and are present in about 20% of cases. ⁶ These mutations occur at conserved arginine residues within the enzyme’s active site, specifically at R132 for IDH1 and at R140 or R172 for IDH2. ⁷ Mutant IDH1 and IDH2 enzymes alter the Krebs cycle and convert isocitrate to 2-hydroxyglutarate (2-HG) instead of alpha-ketoglutarate, leading to the accumulation of 2-HG, which acts as an oncometabolite and disrupts cellular differentiation and proliferation. ⁸ However, there are no reported cancer-associated mutations in IDH3. ⁹ 2-HG inhibits α-KG-dependent enzymes like TET and lysine demethylases, acting as an oncogenic metabolite and causing a hyper-methylated phenotype. ¹⁰ Its accumulation leads to abnormal histone and DNA methylation, blocks cell differentiation, and promotes leukemogenesis. ¹¹ Thus, targeting mutant IDH enzymes is crucial for AML therapy. The development of IDH inhibitors has been a significant advancement in AML therapy. The US Food and Drug Administration (FDA) has approved several IDH inhibitors based on clinical trials demonstrating promising responses. However, the efficacy of these inhibitors compared with standard treatments remains under debate. ¹² Consequently, there is an urgent need for individualized treatments with higher efficacy and better safety profiles.

Mutant IDH1 and IDH2 share ~69% sequence identity, particularly in conserved allosteric pockets at the dimer interface, enabling dual inhibition by compounds like vorasidenib (AG-881) without targeting sequence-unique catalytic residues. ¹³ These pockets, present in both wild-type and mutant forms, form a hydrophobic cleft involving conserved residues such as Val255 in IDH1 and Val294 in IDH2. ¹⁴ Structural studies (eg, IDH1-R132H/NADPH/AG-881 and IDH2-R140Q/NADPH/AG-881 complexes) show that AG-881 binds with nanomolar affinity to the same allosteric site in both isoforms via conserved hydrophilic and hydrophobic interactions, including stacking interactions. Binding induces conformational changes, like α10 helix displacement, that sterically block substrate binding and suppress D-2-HG production by locking the enzyme in an inactive open form. Unlike isoform-selective inhibitors (eg, ivosidenib or enasidenib), which exploit pocket shape differences and regulatory segment stability, AG-881 targets conserved structural motifs and topography. ¹⁵ Subtle isoform differences affect kinetics but do not impair binding. Supported by trials (NCT02481154, NCT02492737) and structural, biochemical, and clinical data, dual inhibitors exploit conserved pocket architectures, not catalytic-site chemistry. These findings are reinforced by global research efforts and funding (eg, NSFC, China’s 973 Program), highlighting the clinical relevance of pan-IDH inhibition.

Natural products remain a prolific source of therapeutics: their stereochemical complexity, dense 3-D architectures and “privileged” scaffolds sample chemical space often inaccessible to purely synthetic libraries, which can translate into high target affinity and favorable safety—factors that have yielded many anticancer drugs.^16,17 To leverage these advantages, we screened 407 270 COCONUT bioactives for dual inhibition of mutant IDH1-R132H and IDH2-R140Q—an important unmet need because clinically advanced dual IDH inhibitors to date (eg, vorasidenib/AG-881) are synthetic and well-characterized dual inhibitors of natural origin remain scarce. We therefore applied an integrated in silico pipeline (ligand-centric pharmacophores, structure-based virtual screening, docking, MM-GBSA), followed by ADMET filtering and 500-ns MD to prioritize stable, drug-like dual-inhibitor candidates. ¹⁸

Our integrated approach aims to identify potent dual inhibitors of IDH1 and IDH2 mutations, contributing to the development of more effective treatments for AML. By leveraging natural products and computational techniques, this research provides valuable insights into the discovery of novel therapeutic agents for leukemia and other IDH-related malignancies. We have graphically described the methodology used in Figure 1.

Figure 1.

The schematic diagram of the complete data set screening steps used in this study.

Materials and Methods

Data Set generation and screening

The COCONUT (version 2022, Collection of Open Natural ProdUcTs Online) database was utilized as the screening source, containing approximately 407 270 molecules. ¹⁹ The compounds underwent refinement through the QikProp tool within the “Structure Screening” setting of Schrödinger 2023-1, ensuring absolute adherence to Lipinski’s Rule of Five, yielding a final count of 276 409 entities. Subsequently, these screened molecules were processed utilizing the LigPrep function, employing the OPLS4 force field. ²⁰ The potential molecular states were produced through Epik [49] within a pH range of 7.0 ± 2.0, preserving designated stereochemistry. This was subsequently followed by the creation of tautomers and the derivation of a maximum of 16 energetically favorable conformers for each ligand. The cumulative count of resulting structures was approximately 3.57 million.

Pharmacophore modeling/screening

Pharmacophore hypothesis generation

A pharmacophore consists of unique traits such as spatial constraints, hydrogen bond donor (HBD), hydrogen bond acceptor (HBA), and electronic molecular properties. These elements define the specific manner in which a molecule interacts within a biological target site. ²¹ The construction of ligand-based pharmacophore models was based on recognized active compounds with validated therapeutic effects against the selected receptor.

For pharmacophore model generation, a preliminary selection of the top 100 compounds associated with IDH1 and IDH2 was retrieved from the binding database and arranged according to their IC₅₀ values. Three approved drugs that are effective in IDH1 inhibition, namely, Enasidenib (89683805), Olutasidenib (118955396), and Ivosidenib (71657455) was employed alongside the corresponding conformers for reference data sets for IDH1 and Enasidenib (89683805) for IDH2. ²² This meticulously assembled data set served as the basis for later pharmacophore modeling and selection. The LigPrep tool utilized the OPLS4 force field during processing, ²⁰ aligned with standard settings for ionization, additional automated steps, including the elimination of salts, stabilization of tautomeric forms, and computational refinements, were executed as per built-in configurations. Well-optimized three-dimensional conformations of druggable compounds were generated.

Within the Develop Pharmacophore Model module, the similarity criterion for the conceptual framework was adjusted to three-fourths in both cases, while the count of essential characteristics within the model ranged between 4 and 7, with an optimal selection of 5. The evaluation and prioritization of the generated models adhered to the predefined “Phase Hypo Score.” ²³ The Generate Conformer and Minimize Output settings were enabled, specifying a target of 50 conformers. ²⁴

Validation of the pharmacophore model

The constructed pharmacophore framework underwent validation to assess its proficiency in accurately forecasting the potency of novel molecules. This step is essential prior to utilizing the model for computational screening. ²⁵ The parameters employed to evaluate the efficiency of the developed pharmacophore model include the enrichment factor (EF), receiver operating characteristic (ROC) curves, Boltzmann-enhanced discrimination of ROC (BEDROC), and robust initial enhancement (RIE). ²⁶ These metrics were selected because ROC/AUC provides overall discrimination, EF gauges how well actives appear early in ranked lists, and BEDROC/RIE apply exponential weighting to further emphasize those top ranks, which is critical in virtual screening where early hits save time and resources.

A collection of misleading compounds was assembled utilizing the DUDE Decoy Generation tool, accessible via the link—dude.docking.org/generate. ²⁷ Open Babel v2.4.1 was employed to transform the output into three-dimensional structures. ²⁸ Ligand processing was carried out using the standard parameters and established procedures of LigPrep. The OPLS4 force field ²⁰ was used for the minimization.

Screening using pharmacophore model

The processed ligands underwent screening with the top-performing pharmacophore models, DDHR_1 for IDH1 and ADDHR_1 for IDH2, utilizing the “Phase Module.” The option to prioritize partial matches with a higher number of features was enabled. All pharmacophore characteristics (4 out of 4 for IDH1 and 5 out of 5 for IDH2) were considered during screening. This process filtered approximately 3.57 million structures down to ~181 K for IDH1 and ~67 K for IDH2.

Molecular interactions

Protein preparation

The X-ray crystallographic structures of the IDH1 and IDH2 target proteins, designated by PDB IDs 5DE1 and 5I96, were sourced from https://www.rcsb.org, featuring resolutions of 2.25Å and 1.55Å, respectively, and utilized in this study. ²⁹ Structural analysis was performed utilizing the PDBsum web server (https://www.ebi.ac.uk/thornton-srv/databases/pdbsum). ³⁰ The PDBsum online server was additionally employed to validate IDH1/2 structures using the Ramachandran plot.

Within the Maestro suite by Schrödinger, the protein refinement tool facilitated initial processing of the proteins utilizing the PROPKA algorithm to enhance hydrogen bond configurations. ³¹ Subsequently, the geometries were optimized until the heavy atoms attained convergence at an RMSD of 0.3Å, employing the OPLS4 potential field, and later, solvent molecules situated beyond a 5Å radius from the ligands were excluded. ²⁰

Receptor grid generation

The receptor grid was generated keeping the hydrophobic region and the region where co-crystallized inhibitor GSK321A and AG-221 (Enasidenib) for IDH1 and IDH2, respectively, was attached to the complex, at the centroid of the grid. The coordinates of the receptor grid for IDH1 were X = −21.12, Y = 11.78, Z = 6.25 and for IDH2 were X = −1.53, Y = 15.61, Z = −24.69.

Molecular docking

Molecular docking was constrained to compounds containing a maximum of 100 flexible torsions and an atomic count below 500. The scaling coefficient for Van der Waals interactions was adjusted to 0.80, while the threshold for partial charges was fixed at 0.15. ³² Nitrogen inversion instances and cyclic structure variations were initiated, while ligand flexibility was incorporated into the sampling process. Every anticipated functional moiety underwent torsional bias adjustments. The system was optimized to enhance internal hydrogen bonding and maintain the coplanarity of conjugated π-systems.

The molecular docking process was carried out using a sequence of progressively stringent filters, starting with the HTVS approach (rapid virtual screening) to swiftly evaluate extensive chemical libraries comprising millions of molecules. This was followed by the SP method (moderate precision), which accurately docked tens to hundreds of thousands of ligands. Finally, the XP strategy (enhanced accuracy) was employed to refine the selection by eliminating potential false positives through comprehensive conformational sampling and sophisticated scoring, leading to superior enrichment. At each stage, only the highest-ranked 10% of candidates advanced to the next level. ³³ The HTVS filtered the library to 1 40 350, SP filtered to 25 872 and, XP filtered to around 2589 structures for IDH1 and 4865, 475, and 49, respectively, for IDH2.^{34 -37}

MM-GBSA screening

The prime module was employed to estimate the energetic values derived through the MM-GBSA (Molecular Mechanics–Generalized Born Surface Area) approach. This analysis was conducted to determine the stabilization energy contributions resulting from the possible interactions between the chosen compounds and the receptor proteins 5DE1 and 5I96. ³⁸ The variable-dielectric generalized born (VSGB) solvation model was used and the force field was set as OPLS4. ³⁹

The highest-ranking 10% of ligands obtained from XP docking were evaluated using MM-GBSA for IDH1, while all ligands generated through XP docking were analyzed for IDH2. And again, the compound having binding affinity with both the proteins was selected from the top 10% of MM-GBSA results. Only the top ~10% of XP-docked ligands were forwarded to MM-GBSA rescoring to balance computational expense with retention of high-confidence hits, a practice widely adopted in virtual screening workflows. The only compound that matched the criteria was CNP0166496 (Ternstroside D) which was used for further ADME-Tox and MD Simulation studies.

In silico ADME/T and toxicity analysis

QikProp Module of Schrödinger, SwissADME, ⁴⁰ ProTox-3.0, ⁴¹ ADMETlab 2.0, ⁴² and Deep-PK ⁴³ servers were utilized for pharmacokinetic analysis to evaluate the comprehensive ADMET characteristics of CNP0166496 (Ternstroside D).

Molecular dynamics (MD) simulation studies

The MD simulations for the 5DE1-CNP0166496 and 5I96-CNP0166496 complexes were executed utilizing the Desmond software suite. Every complex was separately embedded within a 10 Å orthorhombic aqueous environment, employing the TIP3P solvation model. ⁴⁴ The ligand-protein complexes were modeled by the OPLS4 force field. ²⁰ Sodium counterions were incorporated into the protein-ligand complexes to achieve charge neutrality within the systems designated for MD simulations. In addition, energy minimization was performed for 2000 iterations to attain the lowest possible energy state before commencing the simulation under an NPT ensemble pathway. ⁴⁵ To enhance statistical reliability, each simulation was run for 500 nanoseconds in triplicate, yielding 3 independent trajectories per complex. Triplicate 500-ns MD runs were performed to ensure reproducibility across stochastic variables (eg, initial velocities), increasing confidence in stability and binding mode assessments. The average values from triplicates were used in figure generation, while the standard deviation-based error bars were included in the supplementary figures.

The root-mean-square deviation (RMSD) predominantly indicates how stable the ligand remains during interaction, whereas the root-mean-square fluctuation (RMSF) characterizes the dynamic variations and adaptability of amino acid residues in the protein, especially in the key active region essential for pharmaceutical research. Root-mean-square deviation is calculated by the discovery $RMS{D}_x=\sqrt{\frac{1}{N}\sum _{i=1}^N\left(r_i^{\prime }\left(t_x\right)\right)-r_j{\left(t_{ref}\right)}^2}$ . Root-mean-square fluctuation values contribute to characterizing protein structures by offering insights into localized conformational variations within protein chains. Expressed in Ångströms (Å), RMSF is determined using the following equation: $RMS{F}_i=\sqrt{\frac{1}{T}\sum _{t=1}^T<\left(r_i^{\prime }\left(t\right)\right)-r_j\left(t_{ref}\right){)}^2>}$ ⁴⁶

Conformational dynamics and energy landscape analysis

The path data produced following the MD computation was transformed into .dcd structure using VMD 1.9.3 application, ⁴⁷ which was then submitted to the MDM-TASK web server for executing principal component analysis (PCA) and dynamical cross-correlation matrix (DCCM) analyses. ⁴⁸ Principal component analysis was performed to detect significant structural shifts and key movements within the protein-ligand assembly, whereas dynamic cross-correlation mapping evaluated the synchronized actions of amino acid units across the simulation duration. The free energy landscape (FEL) was constructed based on the first 2 principal components using the equation:

$$F(X)=-k_{B}TlnP(X)$$

where F(X) represents the free energy, k_B is the Boltzmann constant, T is the absolute temperature, and P(X) is the probability distribution of the system along the principal components. ⁴⁹

Binding free energy analysis

To assess ligand binding strength, MM-GBSA analysis was conducted on snapshots extracted from the 500-ns triplicate MD trajectories. For choosing protein-ligand pairings, the affinity values derived from this technique proved more effective than glide score metrics. Primary energy components, including hydrogen bond interaction strength (ΔG_{Bind_Hbond}), solvation electrostatic free energy (ΔG_{Bind_Solv}), Coulombic or charge-based interaction energy (ΔG_{Bind_Coul}), hydrophobic interaction energy (ΔG_{Bind_Lipo}), and van der Waals force energy (ΔG_{Bind_vdW}), were collectively factored into determining the MM-GBSA-derived relative binding strength.

Results

The IDH proteins

Evaluation using the PDBsum web server revealed 18 α-helices, 3 β-hairpins, 2 β-bulges, 17 β-turns, and 20 helix-helix interactions for 5DE1 (IDH1), while 5I96 (IDH2) exhibited 18 α-helices, 3 β-hairpins, 1 β-bulge, 18 β-turns, and 25 helix-helix interactions.

In addition, protein validation was conducted using the Ramachandran plot, which revealed that 92% of residues were located in favored regions, 7.4% in allowed regions, 0.6% in generously allowed regions, and none in disallowed regions for IDH1. For IDH2, the distribution was 91.7%, 7.3%, 1.0%, and 0.0%, respectively. The computed overall G-Factor was 0.26 for IDH1 and 0.15 for IDH2. ⁵⁰ (Figure 2).

Figure 2.

The folding pattern of the amino acid chain in protein 5DE (produced via PDBsum); B: Ramachandran diagram for the protein, where red highlights low-energy zones, yellow marks permitted areas, pale yellow denotes broadly acceptable regions, and white signifies restricted zones; C: The three-dimensional form of the IDH1 (Chain B) protein, showing interaction with the co-crystallized compound GSK321A; D: The structural arrangement of the amino acid sequence in protein 5I96 (created using PDBsum); B: Ramachandran chart of the protein; C: The spatial configuration of the IDH2 protein, illustrating engagement with the co-crystallized molecule AG-221 (Enasidenib).

Pharmacophore modeling/screening

The constructed pharmacophore frameworks were systematically ordered according to their survival rates and BEDROC values. For IDH1, the model DHRR_1 exhibited the top Phase Hypo Score of 1.328, featuring 1 hydrogen bond donor (D), 1 hydrophobic region (H), and 2 aromatic ring elements (R). Its survival value reached 6.078, with a site value of 0.779, vector value of 0.933, volume value of 0.634, selectivity value of 1.449, and BEDROC value of 0.963 (Figure 3A and (B). For IDH2, the model ADDHR_1 achieved the highest Phase Hypo Score of 1.330, incorporating one hydrogen bond acceptor (A), 2 hydrogen bond donors (D), one hydrophobic area (H), and one aromatic ring component (R). Its survival value was 6.248, with a site value of 0.837, vector value of 0.943, volume value of 0.723, selectivity value of 1.703, and BEDROC value of 0.956 (Figure 3C and (D).

Figure 3.

The distance (A) and angles (B) between pharmacophore features of the best pharmacophore model, DHRR_1 for IDH1 and ADDHR_1 (C and D) for IDH2.

Approximately 3.5 million processed ligand structures were filtered using the DHRR_1 pharmacophore model, reducing the set to 181 118 compounds for IDH1, while the ADDHR_1 model refined the selection to 67 859 structures for IDH2.

Pharmacophore models validation

The constructed pharmacophore framework was rigorously assessed to evaluate its accuracy in predicting the potency of novel compounds identified through database screening or newly designed molecules. Validation serves as a crucial step before utilizing a pharmacophore model for computational screening. ²⁵

The DHRR_1 model was evaluated using a data set of 200 inactive compounds, which expanded to 551 structures following ligand preparation via the DUDE Decoys tool. This model effectively reduced the decoy set to 35 structures, achieving an efficiency rate of 93.648%. Further validation demonstrated an enrichment factor (EF) of 15.93% within the top 1% of the data set, signifying a 15.93-fold enhancement in identifying active compounds across the entire collection. The ROC score, RIE, and the area under accumulation curve (AUAC) values were calculated as 0.77, 8, and 0.82, respectively, confirming the model’s strong capability to distinguish active compounds from decoys with statistical reliability. To further validate its predictive strength, BEDROC values were computed, emphasizing early active compound detection. BEDROC scores at tuning parameters α = 8.0, α = 20.0, and α = 160.9 were determined as 0.686, 0.702, and 0.995, respectively, reinforcing the model’s statistical robustness (Supplementary Figure 1A).

Likewise, the ADDHR_1 model was applied to the same data set, refining it to 37 structures, demonstrating an efficiency rate of 93.285%. Validation of the ADDHR_1 framework revealed an enrichment factor (EF) of 25.39% within the top 1% of the decoy data set, signifying that the model was 25.39 times more effective in identifying true positives or active compounds from the entire data set. The ROC score, RIE, and AUAC values were computed as 0.94, 12.04, and 0.94, respectively, confirming the statistical reliability of ADDHR_1 in distinguishing active compounds from decoys. To further validate its predictive strength, BEDROC values were calculated, emphasizing early active compound detection. BEDROC scores at tuning parameters α = 8.0, α = 20.0, and α = 160.9 were determined as 0.900, 0.870, and 0.983, respectively, reinforcing the model’s statistical significance (Supplementary Figure 1B).

Binding affinity predictions

The docking scores for the top 500 hit compounds ranged from −16.131 to −11.628 kcal/mol for IDH1, after XP docking, and −16.822 to −7.017 kcal/mol for the 49 compounds for IDH2. The compound having binding affinity with both the proteins was in the top MMGBSA results was CNP0166496 (Ternstroside D) with ΔG_bind energy of 84.45 kcal/mol and −60.73 kcal/mol for IDH1 and IDH2, respectively. It had a molecular formula of C₂₂H₂₆O₁₁.

Ligand CNP0166496 exhibited a variety of interactions with multiple amino acid residues within the binding pockets of both IDH1 and IDH2. It showed major H-bond interactions with Leu120, Lys126, and Ile128 at a distance of 2.20 Å, 1.70 Å, and 1.99 Å, respectively, with IDH1 and Gln(B)316, Asp(B)312, Gln(A)316, Val(A)294, and Asp(A)312 with IDH2, at a distance of 1.85 Å, 2.10 Å, 2.04 Å, 1.95 Å, and 1.99 Å, respectively (Table 1).

Table 1.

Top interacting residues of CNP0166496 (Ternstroside D) with the IDH proteins, XP docking score (kcal/mol) and MM-GBSA ΔG_bind energy (kcal/mol).

Protein	H-bond Interacting Residues	Docking Score	MM-GBSA (ΔGbind)
IDH1	L120, K126, I128	−14.200 kcal/mol	−84.45 kcal/mol
IDH2	Q(B)316, D(B)312, Q(A)316, V(A)294, D(A)312	−16.822 kcal/mol	−60.73 kcal/mol

ADME/T analysis

CNP0166496 shows significant pharmacokinetic and pharmacological potential as an excellent drug candidate. Its high melting and boiling points indicate stability, and its hydration free energy suggests good solubility. With pKa values ensuring stability across pH levels, a dipole moment reflecting moderate polarity, and a topological polar surface area (TPSA) indicating strong biological interactions, it balances structural complexity and adaptability. The compound’s lipophilicity supports effective membrane permeability, and compliance with Lipinski’s Rule of 5 underscores its drug-likeness. Despite minor synthetic and reactivity alerts, it remains viable. Moderate water solubility and Caco-2 permeability are offset by reasonable human intestinal and oral absorption. As a P-glycoprotein substrate, it shows extensive tissue distribution and substantial plasma protein binding. Low blood brain barrier (BBB) and central nervous system (CNS) permeability reduce CNS side effects. Lacking major cytochrome P450 (CYP) enzyme interactions minimizes drug-drug interactions, and efficient clearance prevents accumulation. Its satisfactory safety profile includes predicted GHS class 5 toxicity, moderate AMES toxicity, no hERG inhibition, moderate oral rat toxicity, and no skin sensitization, drug induced Liver Injury (DILI), carcinogenicity, or eye and respiratory toxicity. A low FDA recommended daily dose suggests high potency. Overall, the compound’s stability, favorable pharmacokinetics, and safety profile make it an excellent drug candidate. The physiochemical, medicinal chemistry, and ADMET properties of the compound is documented in Table S1, S2, and S3, respectively. The bioavailability radar plot of CNP0166496 (Ternstroside D) obtained from SwissADME is shown in Figure 4.

Figure 4.

The bioavailability radar plot, generated using SwissADME, illustrates the strong drug-likeness of CNP0166496 (Ternstroside D). The reddish-brown region highlights the optimal range for each pharmacokinetic property. A more comprehensive radar plot, obtained from ADMETlab, is provided in Supplementary Figure S1.

MD simulation

RMSD of IDH proteins

The RMSD analysis of IDH1 and IDH2 proteins bound to CNP0166496 over 500 ns triplicate MD simulations reveals distinct stability profiles (Figure 5, Panel A i and B i). For IDH1, the RMSD plot shows an initial sharp rise within the first 50 ns, reaching ~5.5 Å, indicative of equilibration as the protein adjusts to the ligand-bound state. This is followed by a gradual increase and stabilization around 6.0 ± 0.3 Å, with consistent trends across all replicates. Although the elevated RMSD reflects moderate structural deviation from the starting conformation, the convergence and narrow post-equilibration variance suggest the system reaches a semi-stable state. This stability, despite the higher RMSD, indicates the global fold is preserved while accommodating ligand-induced conformational flexibility, potentially linked to allosteric modulation or domain movement in IDH1. In contrast, the IDH2 RMSD profile exhibits a more compact and stable trajectory, with a brief equilibration phase (25–30 ns) followed by stabilization within 2.0–2.5 ± 0.2 Å. The low and consistent RMSD across all replicates underscores high structural integrity and reproducibility, suggesting that CNP0166496 does not disrupt IDH2’s native conformation. This reinforces the ligand’s compatibility with the IDH2 binding pocket and supports the use of this complex for downstream analyses such as binding free energy calculations or per-residue energetic decomposition.

Figure 5.

Molecular dynamics simulation analysis of IDH1–CNP0166496 (A) and IDH2–CNP0166496 (B) complexes over 500 ns. Panels i and ii represent the RMSD of the protein backbone and ligand, respectively. Panels iii and iv show the RMSF of protein residues and ligand atoms. Panel v displays the number of hydrogen bonds formed over time. Panel vi presents the radius of gyration, while panels vii and viii depict the solvent accessible surface area and polar surface area, respectively.

RMSD of CNP0166496

In the IDH1 complex, the ligand RMSD rises from ~0.5 to ~2.2 Å during the first 100 ns, reflecting early-stage conformational adaptation within the binding site (Figure 5, Panel A ii). Beyond this phase, the RMSD fluctuates moderately between 1.8 and 2.8 Å, with occasional spikes reaching ~3.2–3.5 Å, indicating a semi-flexible ligand conformation marked by dynamic reorientation of certain moieties. Despite this mobility, the average RMSD stabilizes at 2.3 ± 0.4 Å, and the absence of significant drift or detachment suggests the ligand remains consistently bound. Replicate convergence reinforces the robustness of this behavior. In contrast, the ligand in the IDH2 complex demonstrates markedly lower mobility, with initial fluctuations subsiding within 50 ns, followed by stable RMSD values in the narrow range of 0.9–1.2 ± 0.1 Å throughout the simulation (Figure 5, Panel B ii). This low-amplitude variation and tight overlap across replicates suggest a highly stable binding mode, with minimal conformational reorganization or positional drift.

Protein RMSF

For IDH1, the RMSF profile shows that most residues fluctuate within 1.0–2.5 Å, indicating stable secondary structures; however, prominent peaks are observed around residues 125–175 and 250, with fluctuations reaching up to ~9.0 Å. These high-RMSF regions likely correspond to solvent-exposed loops or coil structures that are not directly involved in ligand binding, suggesting surface-level flexibility. Importantly, residues comprising the ligand-binding cavity show minimal fluctuations, implying stable interactions with the ligand. The close alignment across all triplicates confirms the reproducibility of these dynamic patterns. In contrast, the IDH2 RMSF profile is more restrained, with most residue fluctuations confined to 0.8–2.5 Å. Mild peaks at ~75–100, ~150, and ~325 reach just over 5.0 Å and are likely localized to peripheral or loop regions. The core structure remains compact and stable, with ligand-interacting residues exhibiting minimal movement. This rigidity supports the observed lower ligand RMSD and suggests that the IDH2 binding pocket provides a structurally stable and less dynamic environment. The strong agreement among replicates further confirms the reliability of this system (Figure 5, Panel A iii and B iii). All the corresponding averaged data with standard deviation error bars is provided in Supplementary Figure S3 and Supplementary Figure S4.

Ligand RMSF

In the IDH1–CNP0166496 complex, the RMSF values gradually increase from ~1.0 Å at the central atoms to a peak of ~2.5 Å at atoms 30–33, indicating enhanced flexibility at the ligand’s termini—likely due to solvent exposure or reduced steric constraint within the binding cavity. The central region (atoms ~5–15) exhibits lower RMSF values of 0.8–1.5 Å, suggesting a well-anchored core scaffold. The modest variability across replicates affirms reproducible dynamics and supports the interpretation that while the ligand’s core remains stably engaged, peripheral regions undergo limited movement, allowing conformational plasticity potentially important for adaptive binding in IDH1. In contrast, the IDH2–CNP0166496 complex exhibits a more rigid ligand profile, with RMSF values consistently between 0.5 and 0.75 Å and only a minor peak at atom 20 (~1.3 Å). This tightly constrained fluctuation pattern indicates a stable, well-locked ligand conformation across all atoms, reflecting a structurally restrictive and complementary binding pocket in IDH2. The narrow overlap across triplicates further confirms reproducibility and minimal conformational drift (Figure 5, Panel A iv and B iv).

Binding dynamics of CNP0166496 in IDH1 and IDH2

The hydrogen bond analysis of CNP0166496 bound to IDH1 and IDH2 reveals contrasting yet functionally complementary interaction landscapes, shedding light on their respective binding dynamics over 500 ns MD simulations (Figure 5, Panel A v and B v). In the IDH1–CNP0166496 complex, a dynamic yet persistent hydrogen bonding network is observed, with the ligand forming an average of 2–3 hydrogen bonds throughout the simulation and occasionally rising to 4–5, particularly between 250 and 450 ns. These interactions—especially with residues such as VAL121, LEU120, ILE128, and the catalytically relevant TYR285—occur across the polyhydroxylated core and terminal hydroxyl groups of the ligand, often involving both direct and water-mediated bridges. The contact fraction plot highlights VAL121 as the most interactive residue, with additional stabilizing contributions from ILE126 and GLU132, suggesting multivalent anchoring that allows the ligand to remain stably engaged while dynamically sampling different microenvironments within the binding pocket. This flexible yet stable binding mode indicates strong affinity and specificity, driven by transient hydrogen bonding hotspots (Figure 6, Panel A). In contrast, the IDH2–CNP0166496 complex exhibits a more spatially coordinated and less fluctuating hydrogen bond pattern. Although the average hydrogen bond count remains low (~0.3), key polar interactions with residues ASP312, GLN316, and VAL294 provide critical anchoring, particularly with ASP312 forming persistent dual-chain interactions (>90% occupancy) that span both subunits—a unique cross-subunit stabilization not seen in IDH1. Contact dynamics and fraction plots confirm sustained engagement with GLN316, VAL315, and ILE319, and additional hydrophobic support from VAL294 and VAL297 further reinforces a deeply embedded, enthalpically favorable binding pose. The low ligand RMSD and RMSF values in IDH2 support this interpretation, suggesting that despite fewer hydrogen bonds, the ligand experiences enhanced shape complementarity and entropic stabilization, contributing to tight, persistent binding (Figure 6, Panel B).

Figure 6.

Protein–ligand interaction analysis of IDH1–CNP0166496 (A) and IDH2–CNP0166496 (B) complexes. Panel A i and B i depict the protein–ligand contact dynamics showing key interacting residues with the ligand CNP0166496 over the simulation trajectory. Panel A ii and B ii illustrate the interaction fraction profiles, indicating the frequency and types of interactions (hydrogen bonds, hydrophobic, ionic, and water bridges) between the ligand and protein residues throughout the 500 ns molecular dynamics simulation.

Radius of gyration (Rg)

In the IDH1 complex, the Rg fluctuates moderately between ~4.2 and 5.2 Å, with an average around 4.6 ± 0.2 Å. A gradual decrease in Rg values over time suggests progressive structural compaction, while intermittent spikes—such as those around 200 ns (~5.5 Å) and 400 ns (~5.3 Å)—point to transient expansion events, possibly driven by loop flexibility or partial relaxation of peripheral regions (Figure 5, Panel A vi). Nonetheless, the return to lower Rg values and tight triplicate convergence indicates that the complex regains compactness and maintains overall conformational stability. This dynamic yet convergent profile suggests subtle reorganization favoring a tighter, more stable fold. In contrast, the Rg profile of the IDH2–CNP0166496 complex remains exceptionally stable throughout the trajectory, with values tightly centered around 4.95 ± 0.05 Å and minimal variation across replicates. No major expansion or contraction events are observed, underscoring the rigid, well-conserved structure of the complex. This structural integrity aligns with previously noted low protein RMSD and ligand RMSF values, reinforcing the conclusion that IDH2 forms a compact, reproducibly stable assembly with CNP0166496 (Figure 5, Panel B vi).

Solvent accessibility and polar surface analysis

In the IDH1 complex, solvent accessible surface area (SASA) values increase gradually during the first 200 ns—from ~100 Ų to peaks of ~200–250 Ų, with occasional spikes exceeding 300 Ų—indicating transient surface expansions, possibly from loop motions or partial ligand exposure. Beyond 300 ns, the SASA stabilizes between 120 and 150 Ų, suggesting the system settles into a more compact, less solvent-exposed conformation. These fluctuations, while pronounced early, do not compromise overall complex stability and likely reflect localized rearrangements near the solvent interface. By comparison, the IDH2 complex exhibits a more restrained SASA trajectory, with values plateauing around 150–200 Ų after 200 ns and remaining within a steady envelope of 130–160 Ų. Although occasional spikes (up to 400 Ų) occur, such as around 425 ns, the overall solvent exposure is more consistently controlled than in IDH1, aligning with IDH2’s lower RMSD and ligand flexibility (Figure 5, Panel A vii and Panel B vii).

In terms of PSA, IDH1 exhibits moderate fluctuation within a consistent range of 320–360 Ų (average ~340 ± 10 Ų), suggesting sustained polar interactions and solubility potential. Occasional spikes (~370 Ų) and dips (~310 Ų) indicate transient conformational adjustments or brief exposure of polar groups, but the tight clustering across replicates confirms robustness and reproducibility. This supports the hypothesis that the ligand retains hydrogen-bonding capacity and a flexible, aqueous-compatible profile in the IDH1 pocket. In contrast, the IDH2 complex shows a narrower, more elevated PSA band, fluctuating around ~368 ± 5 Ų with minimal deviation (360–375 Ų), reflecting exceptional polar surface stability and consistent solvent accessibility of the ligand’s polar moieties. Despite fewer direct H-bonds with IDH2, the high PSA suggests dominant solvation-based polarity, likely enhancing bioavailability (Figure 5, Panel A viii and Panel B viii).

Conformational dynamics

The PCA and FEL analyses of the IDH1– and IDH2–CNP0166496 complexes provide contrasting insights into their conformational dynamics over 500 ns MD simulations. For the IDH1 complex, PCA reveals 2 dominant clusters in PC1–PC2 space, indicating transitions between 2 major energy basins. The pronounced variance along PC1 (>10 units) highlights significant conformational shifts, likely stemming from loop or domain-level motions. Correspondingly, the 3D FEL plotted against Radius of Gyration (Rg) and RMSD displays 2 well-defined, energetically favorable minima, suggesting the presence of metastable conformational states sampled by the complex. The 2D FEL further reinforces this with a distinct basin localized around RMSD ~2.5–3.5 and Rg ~21–24 Å, indicating that the complex stabilizes within a specific conformational zone, undergoing localized but not global structural transitions. This suggests that IDH1, when bound to CNP0166496, maintains consistent structural integrity while transitioning between limited, stable conformers (Figure 7, Panel A).

Figure 7.

PCA and FEL analysis of IDH1–CNP0166496 (A) and IDH2–CNP0166496 (B) complexes during 500 ns molecular dynamics simulation. Panels A i and B i represent the PCA projections along the first 2 principal components (PC1 and PC2), reflecting the conformational clustering and sampling. Panels A ii and B ii show the 3D Free Energy Landscapes plotted as a function of RMSD and radius of gyration, highlighting energy minima corresponding to stable conformational states. Panels A iii and B iii present the 2D Free Energy Landscapes (FEL) further emphasizing the dominant low-energy basins sampled by the complexes over the trajectory.

In contrast, the IDH2–CNP0166496 complex (Panel B) exhibits a broader and more dispersed PCA distribution (B i) across both positive and negative PC1 values, lacking a dominant clustering pattern. This indicates that IDH2 samples a wider array of conformational states, reflecting greater intrinsic flexibility or a more dynamically adaptive binding interaction. The 3D FEL (B ii) displays multiple shallow wells and a flatter energy surface compared with IDH1, supporting the notion of reduced conformational rigidity and increased sampling of low-energy conformers. The 2D FEL (B) depicts a continuous low-energy region extending across Rg ~20–24 and RMSD ~2–3.5 Å, with smooth transitions rather than discrete basins. This suggests that the IDH2 complex resides within a rugged yet interconnected energy landscape, favoring entropic adaptability over rigid stabilization. Together, these results indicate that while CNP0166496 induces discrete and metastable states in IDH1, it promotes dynamic and broader conformational sampling in IDH2—highlighting differential binding-induced flexibility that may influence functional and pharmacological outcomes (Figure 7, Panel B). The Internal PCA, multi-dimensional scaling (MDS), and t-distributed Stochastic Neighbor Embedding (t-SNE) graphs are shown in Supplementary Figure S5.

Dynamical cross-correlation matrix analysis (Supplementary Figure S6) of CNP0166496-bound IDH1 and IDH2 reveals ligand-induced correlated motions—IDH1 shows strong positive and negative correlations around residues 0–100, ~140–170, and 300–400, suggesting coordinated movements stabilizing the structure, while IDH2 displays dispersed correlations and anti-correlations, notably in residues 100–300 and 500–700 with complex upper-right quadrant patterns, indicating ligand-induced conformational changes aligned with PCA results.

Binding free energy analysis

The parameters of free binding energies between the ligands were evaluated by the MM-GBSA simulations, which is the final step for checking the stability levels of the examined systems in the aqueous environment. The average values of calculated MM-GBSA energy parameters are presented in Supplementary Table S4 and Supplementary Table S5 and contain binding free energies, Coulomb energy (Coulomb), Covalent bonding (Covalent), Hydrogen bonding (H-bond), lipophilic bonding (Lipo), π-π stacking interaction (Packing), solvent generalized bonding (SolvGB), and Van der Waals bonding energy (VDW) over the 500 ns time frame. The MM-GBSA binding free energy profile for the IDH1–CNP0166496 complex (Figure 8, Panel A) reveals a gradual destabilization trend over the 500 ns simulation period. Initially, the ΔG_bind values exhibit strong binding affinity around −85 to −90 kcal/mol, but a progressive increase toward −60 kcal/mol is observed, indicating a relative weakening of the ligand’s binding strength. The fluctuation amplitude across triplicates is also noticeable, particularly in the 350–450 ns window, suggesting that the complex undergoes significant conformational rearrangements affecting the interaction stability. These results align with the PCA and FEL findings, where IDH1 displayed distinct conformational transitions between metastable states. The upward shift in binding free energy possibly reflects a movement away from an initial tightly bound conformation toward a more relaxed or less favorable binding pose, possibly due to induced-fit effects or entropic contributions destabilizing ligand-protein contacts over time.

Figure 8.

Binding free energy analysis of the IDH1–CNP0166496 (A) and IDH2–CNP0166496 (B) complexes computed using the MM-GBSA method over 500 ns simulation. The plots show ΔG_bind values for 3 replicates (gray traces) and their average (red trace) across the simulation time, representing the dynamic stability and strength of ligand binding.

The IDH2–CNP0166496 complex (Figure 8, Panel B) demonstrates a consistently strong and stable binding profile across the entire 500 ns trajectory. The ΔG_bind values remain within a narrow range around −75 to −85 kcal/mol with minimal drift and limited variability between the triplicate runs. This reflects a highly stable binding mode where the ligand maintains persistent and favorable interactions with the IDH2 binding pocket. The comparatively low fluctuation amplitude indicates reduced conformational perturbation within the binding interface, consistent with the broader and interconnected FEL landscape seen in the PCA and FEL plots. Despite the greater conformational flexibility observed in IDH2’s PCA, this result suggests that the binding site retains its structural integrity and accommodates CNP0166496 effectively, potentially due to a more adaptable pocket geometry.

Discussion

This in silico study aimed to identify a dual inhibitor targeting IDH1 and IDH2 mutations implicated in AML. Virtual screening identified a lead compound, CNP0166496 (Ternstroside D), exhibiting high docking scores and binding affinities with both IDH1 (−14.2 kcal/mol) and IDH2 (−16.822 kcal/mol) as predicted by MM-GBSA calculations, which are significantly higher than those of the native inhibitors GSK321A for IDH1 (−9.6 kcal/mol) and Enasidenib (AG-221) for IDH2 (−8.9 kcal/mol). ⁵¹ According to the COCONUT database, Ternstroside D was retrieved from various database collections, including NP-MRD, ChEBI, SuperNatural2, UNPD, and NPASS, highlighting its well-documented presence in natural product databases. Ternstroside D is a phytocompound classified as a triterpenoid saponin, characterized by a steroid-like core structure with glycosidic moieties attached, and is isolated from plants in the genus Ternstroemia within the Theaceae family. Notably, Ternstroemia gymnanthera and Ternstroemia japonica, native to Asia and East Asia, respectively, are known sources of this compound. Ternstroside D is a β-D-glucoside featuring a 2-(3,4-dihydroxyphenyl)ethoxy residue at the anomeric position and a [(3,4-dihydroxyphenyl)acetyl]oxy residue at position 6. Notably, its antioxidant activity has been validated in vitro. ⁵² The compound exhibited favorable pharmacokinetic and pharmacological properties, suggesting good drug-likeness according to Lipinski’s Ro5. It demonstrated acceptable water solubility, Caco-2 permeability, and limited blood-brain barrier penetration, indicating potential oral bioavailability and reduced central nervous system side effects. Minor synthetic challenges and potential reactivity issues could be addressed through medicinal chemistry optimization. ⁵³ RMSD analysis indicated stable protein-ligand complex formation throughout the 500 ns simulations for both IDH1 and IDH2. For the IDH1 complex, although the RMSD exceeded the conventional 3 Å threshold, it plateaued after the initial equilibration phase and remained stable throughout the simulation, indicating that the system reached a consistent conformational state. Such elevated RMSD values are common in larger proteins with flexible regions and are widely accepted when the trajectory stabilizes, as observed here.^54,55 While the ligand showed some conformational flexibility, it maintained interactions with the binding pocket. The study identified interactions with residues like Lys126, Tyr285, Trp124, Val281, and Gly284 in IDH1, which are crucial for binding of established IDH1 inhibitors like GSK321A. ⁹ This overlap suggests that Ternstroside D might utilize a similar binding mode, potentially mimicking the inhibitory mechanism of existing drugs. Similarly, interactions with residues like Lys180, Gly145, Try210, His173, and Ala174 in IDH2 were observed, resembling the binding mode of inhibitors like Enasidenib (AG-221), ⁵⁶ suggesting potential competitive binding with the established inhibitor for the same active site. AG-881 (Vorasidenib), a known dual mutant IDH1/2 inhibitor, has been shown via crystallography to bind IDH1-R132H and IDH2-R140Q at a similar allosteric pocket at the dimer interface. In our modeling, AG-881 exhibited docking scores of −8.564 kcal/mol and −13.689 kcal/mol for IDH1 and IDH2, respectively, with corresponding MM-GBSA binding energies of −39.28 kcal/mol and −52.77 kcal/mol. By contrast, Ternstroside D not only engages several of the conserved residues that AG-881 does but also forms additional hydrogen bonds (eg, with Asp312 and Gln316 in IDH2), which is consistent with its superior docking affinities and more favorable MM-GBSA binding energies observed for both IDH1 and IDH2 The presence of hydrogen bonds between Ternstroside D and both IDH1 and IDH2 suggests a crucial role in stabilizing the protein-ligand complex, which could be further analyzed to identify specific functional groups on the ligand contributing to binding affinity. Hydrophobic interactions likely play a role in anchoring the ligand within the active site, with understanding these interactions valuable for future structure-activity relationship (SAR) studies aimed at optimizing potency. The identified interactions with residues such as Lys126, Tyr285, Trp124, Val281, and Gly284 in IDH1 suggest a mechanism similar to established IDH1 inhibitors like GSK321A, while interactions with Lys180, Gly145, Try210, His173, and Ala174 in IDH2 resemble the binding mode of inhibitors like Enasidenib (AG-221), targeting the substrate binding pocket and likely competitively inhibiting the mutant enzyme’s ability to convert isocitrate to 2-HG. While the computational findings are strong, they are subject to inherent limitations of in silico methods. Docking and MM-GBSA approximate solvation and entropic contributions only implicitly, and the force fields used cannot fully recapitulate induced fit or polarization effects of proteins. ⁵⁷ In addition, AG-881 binding data come from allosteric pocket binding confirmed in crystallographic studies, while our predictions for Ternstroside D are yet to be verified in vitro. Future directions might involve evaluating the inhibitory activity of Ternstroside D against cell lines harboring IDH1 and IDH2 mutations to validate the in silico predictions, investigating the mechanism by which Ternstroside D inhibits IDH enzyme activity, assessing the selectivity of Ternstroside D for mutant IDH proteins over other related proteins, and conducting SAR studies to optimize the structure of Ternstroside D to improve binding affinity and potency if possible.

Conclusion

Ternstroside D (CNP0166496), a natural product identified through a comprehensive computational pipeline, emerged as a promising dual inhibitor for both IDH1 and IDH2 mutations. The study’s robust methodology, encompassing pharmacophore modeling and screening, docking, and MD simulations, provides compelling evidence for Ternstroside D’s favorable binding affinity and stable interactions with the target proteins. Furthermore, in silico ADME/T analysis suggests its potential for successful drug development.

These findings pave the way for further exploration of Ternstroside D’s therapeutic potential. In vitro and in vivo studies are crucial next steps to validate the in silico data and translate this discovery into a tangible clinical benefit for AML patients harboring IDH mutations. Moreover, this study underscores the power of in silico methodologies as valuable tools for accelerating drug discovery, particularly in the realm of natural products with immense therapeutic potential.

Supplemental Material

sj-docx-1-bbi-10.1177_11779322251399077 – Supplemental material for Exploration of Natural Products for Targeting IDH1/2 Mutations in Acute Myeloid Leukemia Through Ligand-Based Pharmacophore Screening, Docking, ADME-T, and Molecular Dynamic Simulation Approaches

Supplemental material, sj-docx-1-bbi-10.1177_11779322251399077 for Exploration of Natural Products for Targeting IDH1/2 Mutations in Acute Myeloid Leukemia Through Ligand-Based Pharmacophore Screening, Docking, ADME-T, and Molecular Dynamic Simulation Approaches by Uddalak Das, Dheemanth Reddy Regati, Jitendra Kumar and R Sowdhamini in Bioinformatics and Biology Insights