Integrating machine learning, docking analysis, molecular dynamics, and experimental validation for accelerated discovery of novel FLT3 inhibitors against AML

2.1. Construction of FLT3 inhibitor activity dataset

In this study, a classification model dataset (CHEMBL1974) was sourced from the ChEMBL online database on December 1, 2023 [27]. To ensure accurate classification, the inhibitors collected were divided into active inhibitors (marked as ‘1’) and inactive inhibitors (marked as ‘0’). Generally, IC₅₀ values for approved and clinical FLT3 drugs are below 100 nM. Therefore, a threshold of 100 nM was set in this study: inhibitors with IC₅₀ values under 100 nM were classified as active, while those with values above 1000 nM were deemed inactive. To reduce the influence of boundary effects and experimental errors on predictions, inhibitors with IC₅₀ values between 100 nM and 1000 nM were excluded. For the classification model, we employed 1,189 active inhibitors and 1,412 inactive inhibitors. Using DataWarrior software [28], we analyzed the plain ring systems and Murcko scaffolds of the 2,601 inhibitors, identifying 265 plain rings and 1,141 Murcko scaffolds. Figure 1 highlights the 24 most frequent plain rings, showcasing the compound diversity within the dataset, which contributed to the construction of robust predictive models. Table S1 and S2 provide the SMILES and occurrence frequencies of the top 24 plain rings and Murcko scaffolds after analyzing the 2,601 FLT3 inhibitors.

Close

Figure 1.

The top 24 plain rings of 2601 FLT3 inhibitors.

Export to PPT

The dataset at the cellular level (MV4-11, all compounds are FLT3 inhibitors) for this research was compiled through an extensive literature review and meticulous data collection, resulting in 844 sets of IC₅₀ data for FLT3 inhibitors at the MV4-11 cell level. For compounds with multiple IC₅₀ values, the average value was calculated and used as the final result. In this study, we used the base-10 logarithm of IC₅₀ (expressed as -logIC₅₀, or pIC₅₀) as the dependent variable, rather than IC₅₀. The distribution of this data has been shown in Figure 2, with a notable clustering of pIC₅₀ values between 5.00 and 9.00. DataWarrior software was again employed to examine the plain ring systems and Murcko scaffolds of the 844 FLT3 inhibitors at the MV4-11 cell level, identifying 92 plain rings and 337 Murcko scaffolds. These findings highlight the compound diversity within our MV4-11 dataset. Table S3 and S4 present the SMILES and occurrence frequencies of the top 24 plain rings and Murcko scaffolds after analyzing 844 FLT3 inhibitors at the MV4-11 cell level.

Close

Figure 2.

Distribution of pIC₅₀ data for 844 FLT3 inhibitors at the MV4-11 cell level.

Export to PPT

2.2. ML models building

The molecular fingerprints employed in this study served as inputs for ML algorithms, facilitating the transformation of chemical structures into numerical representations through descriptor computation. These conversions, expressed as either binary or continuous values, enable molecular diversity and similarity to be quantified in a manner suitable for ML models. By encoding molecular data into fixed-length vectors, processing efficiency was substantially improved. Specifically, fingerprints were calculated using the PaDEL software [29], with four types of molecular fingerprints utilized: CDK fingerprint, CDK extended fingerprint, Substructure fingerprint, and Substructure fingerprint count. These fingerprints capture essential molecular properties, aiding in the development of predictive models for FLT3 activity. Additional information on the molecular fingerprints is available in Table S5.

To create binary classification models, LightGBM, Random Forest, k-nearest neighbors (KNN), multi-layer perceptron (MLP), Decision Tree, and Logistic Regression were employed. For regression, LightGBM, Random Forest, KNN, PLS, LASSO, and XGBoost were used. Among these, the LightGBM algorithm was preferred due to its advantages as an ensemble learning method [30]. Additionally, the LightGBM algorithm is capable of building both classification models and powerful regression models. LightGBM, a modification of Gradient Boosting Decision Trees (GBDT), is well-known for its efficiency in tackling classification and regression challenges involving extensive datasets and high-dimensional attributes. It incorporates two primary strategies: gradient-based one-sided sampling (GOSS) and exclusive feature bundling (EFB) [31]. GOSS enables efficient handling of large sample sizes, while EFB optimizes the management of numerous features, reducing overfitting risks. Furthermore, LightGBM’s leaf-wise tree splitting strategy enhances accuracy by focusing on nodes with the greatest delta loss. This approach has consistently outperformed traditional level-wise splitting methods [32]. A ten-fold cross-validation strategy was employed to minimize the impact of random data splitting, and enhancement methods such as hyper-parameter adjustment were utilized to boost the model’s accuracy. The grid search approach was employed to identify the best parameters for the models. When the parameter space is limited and computational resources are sufficient, grid search ensures a comprehensive exploration of parameter combinations and the reproducibility of results, while efficiently utilizing parallel computing to accelerate the optimization process. The datasets were split into training and testing subsets with a 90:10 ratio. In this study, compounds with IC_50≤100 nM were classified as the “positive class,” while those with IC_50≥1000 nM were classified as the “negative class.” The confusion matrix was used to visually represent the model’s performance, comparing true positive (TP), false negative (FN), false positive (FP), and true negative (TN) predictions. Accuracy was defined as (Eq. 1):

When class distribution is balanced, accuracy is a reliable metric, but for imbalanced classes, it tends to favor the majority class. Therefore, the F1-score, which balances precision and recall, was used. Precision and recall were calculated as follows (Eqs. 2-4):

The formulas show that recall depends on the true label, while precision is influenced by predictions. A high recall usually comes with low precision, and vice versa, depending on how the model balances these factors. The F1 score combines both metrics, with higher scores indicating superior model performance.

For regression models, evaluation metrics included Mean absolute error (MAE), f (RMSE), and coefficient of determination (R²). The corresponding formulas are (Eqs. 5-7):

In this equation, y_exp and y_pred represent the true and forecasted values, respectively. ȳ indicates the average of the predicted values, while m denotes the number of samples. The complete source code and dataset are available at https://github.com/Yihuan-Zhao93/FLT3ML.

2.3. Molecular docking studies and molecular dynamics (MD) simulation

To identify novel FLT3 inhibitors, we performed a structure-based virtual screening approach. The crystal structure of FLT3 (PDB ID: 4XUF) was retrieved from the Protein Data Bank (PDB) [33]. Initially, the PDB file was prepared by removing crystallographic water molecules and the bound ligand, followed by assigning partial charges. The 3D structure of the ligand was generated and saved as a pdbqt file using AutoDock software [34]. Once the protein and ligand structures were prepared, molecular docking was conducted with AutoDock Vina. Docking simulations were carried out by placing the ligands into the active site of FLT3, with the grid box set to 20 × 20 × 20 ų, centered at the original ligand’s binding site. Binding affinity was evaluated based on the docking scores, where a more negative binding energy suggested stronger interactions between the ligand and the FLT3 receptor. For detailed interaction analysis, Discovery Studio Visualizer [35] was employed, allowing visualization of the molecular interactions between the ligands and protein. To ensure the reliability of the docking procedure, the co-crystallized ligands were re-docked into their respective binding sites within the FLT3 crystal structure. The resulting docking conformations closely matched the experimentally observed ones, with low root-mean-square deviation (RMSD) values (Figure S1), validating the docking methodology. Subsequently, a virtual screening of approximately 1.6 million compounds from the ChemDiv database was performed using the Glide software. Docking poses and scores were generated for all compounds, providing a comprehensive evaluation of potential inhibitors. Notably, the ChemDiv compounds are commercially available, facilitating rapid procurement and experimental verification.

Utilizing the Gromacs 2022 suite [36], we conducted molecular dynamics (MD) simulations. We followed the detailed protocol as outlined in prior literature [22] to execute a 100 ns molecular dynamics simulation, capturing frames every 10 ps. After concluding the simulation, trajectory analysis was performed utilizing Visual molecular dynamics (VMD) and PyMOL, and the binding free energy between the protein and ligand was calculated using the g_mmpbsa tool.

2.4. Experimental validation

In this study, compounds C1-C4 were procured from TargetMol (Shanghai). Cell viability was assessed using the = (CCK-8) assay. MV4-11 cells were seeded into 96-well plates at a density of 3 × 10⁴ cells per well in 100 μL of medium and incubated for 24 h. The experimental setup included a blank control group, which contained only the 1640 complete medium, and treatment groups with varying concentrations of compounds C1-C4 at 1, 10, 100, 1,000, 10,000, and 100,000 nM. Each treatment condition was replicated in three wells, and the absorbance in each well was recorded three times to ensure accuracy. Compounds C1-C4 were dissolved and diluted in 1640 complete medium to prepare the respective treatment solutions. Following a 24-h treatment period with the compounds, 10 μL of CCK-8 solution was added to each well. Subsequently, the plates were incubated for an additional 2 h at 37°C. Absorbance was then measured at 450 nm using a microplate reader to evaluate cell viability. The percentage of cell viability was calculated with the following formula (Eq. 8):

3.1. Construction of FLT3 biological activity ML classification models

In this study, ML models were developed to predict the activity of FLT3 inhibitors, utilizing various molecular fingerprints as input features. The representation of FLT3 inhibitors in the datasets included CDK fingerprints, CDK extended fingerprints, substructure fingerprints, and substructure fingerprint counts. CDK fingerprints are 1024-bit one-dimensional arrays that capture the presence of distinct structural features. The extended CDK fingerprints improve upon standard CDK fingerprints by incorporating additional ring features. Substructure fingerprints and counts, represented as 307-bit binary strings, indicate the presence of SMiles ARbitrary Target Specification (SMARTS) patterns for functional group classification, as defined by Christian Laggner. Given that model performance can vary significantly based on the choice of input features, we first assessed the predictive performance of models with different input representations. Table 1 presents the performance metrics for various models using ten-fold cross-validation on an independent test set. The results demonstrated that the model employing CDK extended fingerprints (1024 bits) achieved the highest accuracy, precision, and F1-score, indicating superior predictive performance. Combining different molecular fingerprints did not enhance the model’s performance; therefore, CDK extended fingerprints (1024 bits) were selected as the optimal input for the classification models. Figure 3(a) displays the confusion matrix for the classification model using CDK extended fingerprints without cross-validation, with prediction accuracies of 0.959 for active compounds and 0.957 for inactive compounds, resulting in an overall accuracy of 0.958, demonstrating high predictive accuracy. Additionally, the performance of the LightGBM algorithm was compared to five other algorithms. The optimal parameter values for different classification ML prediction models have been shown in Table S6. Figure 3(b) shows the comparative F1-scores of various algorithms utilizing CDK extended fingerprints as input. The LightGBM algorithm exhibited superior predictive performance, validating its effectiveness in constructing structure-activity relationship models for inhibitors. Thus, the LightGBM algorithm was also used in the subsequent activity prediction model at the MV4-11 cell level. To further validate the model’s accuracy, 400 compounds not included in the training or independent test sets (out-of-sample, sourced from the BindingDB database) were selected, consisting of 200 active and 200 inactive compounds. The model accurately predicted the activity of 391 out of the 400 compounds, achieving an overall accuracy of 0.978. The detailed information of the 400 compounds has been listed in Table S7. This high accuracy underscores the model’s robust generalization capability, indicating its potential in predicting FLT3 activity for unknown samples. To illustrate the higher predictive accuracy of our model, we downloaded the data from the referenced study [25] and used CDK extended fingerprints as input, applying the LightGBM method to build a classification model. With ten-fold cross-validation, the training set achieved an accuracy of 0.9782 ± 0.0018, precision of 0.9782 ± 0.0019, and F1 score of 0.9782 ± 0.0018. For the independent test set, the model yielded an accuracy of 0.8427 ± 0.0186, a precision of 0.8437 ± 0.0185, and an F1 score of 0.8426 ± 0.0186. Additionally, we downloaded the 1,350 data points from the recent publication [26] and applied the LightGBM method to build a regression model. With ten-fold cross-validation, the training set achieved an R² of 0.998 ± 0.000, MAE of 0.031 ± 0.001, and RMSE of 0.047 ± 0.001. The test set achieved an R² of 0.896 ± 0.005, MAE of 0.242 ± 0.007, and RMSE of 0.313 ± 0.008. These results further demonstrate the high accuracy of our ML method.

Close

Figure 3.

The performance of FLT3 classification prediction model (a) Confusion matrix (b) Comparison of prediction performance of different algorithms utilizing CDK extended fingerprints as input.

Export to PPT

3.2. Construction of MV4-11 ML models

To identify potent FLT3 inhibitors for targeted AML therapy, we developed ML predictive models based on activity data from 844 FLT3 inhibitors tested in MV4-11 cells. Initially, we evaluated the impact of different molecular fingerprint inputs on prediction performance, as summarized in Table 2. The results show that using individual molecular fingerprints, such as CDK or CDK extended fingerprints, resulted in high prediction accuracy. Specifically, with CDK or CDK extended fingerprints alone, the LightGBM model achieved R² values of 0.809 and 0.805, MAE values of 0.402 and 0.406, and RMSE values of 0.548 and 0.553, respectively, indicating that even single molecular fingerprints provide reasonably accurate predictions. Moreover, combining different molecular fingerprints further improved the model’s performance. For instance, the combination of CDK extended and Substructure count fingerprints yielded an R² value of 0.817, an MAE of 0.394, and an RMSE of 0.536. Notably, when combining CDK, CDK extended, Substructure, and Substructure count fingerprints, the model’s performance reached its highest level, achieving an R² value of 0.819, an MAE of 0.394, and an RMSE of 0.533 on the independent test set. This combination significantly enhanced both prediction accuracy and model stability compared to using single fingerprints alone. Even though integrating all molecular fingerprints yielded the highest predictive performance, utilizing only the CDK, CDK extended, and Substructure count fingerprints resulted in similarly effective outcomes. Consequently, we chose this combination as the ideal input for our ultimate predictive model. In summary, the combination of multiple molecular fingerprints substantially enhances the performance of FLT3 inhibitor activity prediction models. Figure 4 illustrates the prediction performance of the model using CDK, CDK extended, and Substructure count fingerprints as inputs. For the independent test set, the model achieved an R² value of 0.846, an MAE of 0.368, and an RMSE of 0.492, demonstrating excellent predictive capabilities. Figure S2 presents a comparison of the model performance between LightGBM and other algorithms, including Random Forest, KNN, PLS, LASSO, and XGBoost. From the figure, it can be observed that the model constructed using the LightGBM algorithm has the smallest MAE and RMSE. Table S8 presents the parameter values of different ML prediction models. Therefore, after constructing classification and regression prediction models with high prediction accuracy, we applied these models to predict the activity of compounds that exhibited high docking scores obtained through molecular docking.

Close

Figure 4.

The true pIC₅₀ value and predicted value obtained from LightGBM model using CDK, CDK extended, and Substructure count fingerprints as inputs of (a) Training set and (b) Independent test set.

Export to PPT

3.3. Molecular docking studies

To ensure the reliability of molecular docking, we thoroughly validated it beforehand. The ligand P30 (quizartinib) was re-docked into the FLT3 binding site, achieving an RMSD of 1.481 Å, as illustrated in Figure S2, confirming accurate pose prediction. The binding energy between P30 and FLT3 was determined to be -10.524 kcal/mol, indicating a robust interaction. Additionally, we performed docking with sorafenib and FLT3, which yielded a docking score of -10.997 kcal/mol. This supports the use of -10.524 kcal/mol as a reliable threshold for screening potent FLT3 inhibitors. Next, we docked 1,600,000 molecules from the ChemDiv library into the FLT3 site using set parameters. A threshold of -10.524 kcal/mol was established, selecting 7,280 compounds as potential FLT3 inhibitors for further study. This criterion was chosen to identify compounds with strong binding affinity. We employed DataWarrior software to explore the structural diversity, analyzing plain ring systems and Murcko scaffolds. This revealed 436 unique rings and 2,678 scaffolds, indicating notable diversity beneficial for discovering new FLT3 inhibitors. Molecular weight and AlogP values were calculated and plotted, with Figure 5(a) showing a wide chemical space distribution, highlighting significant diversity. ML models predicted compound activity, identifying 2,720 as active, with 304 having a high activity probability, indicating strong potential as FLT3 inhibitors. The MV4-11 prediction model assessed cellular activity, finding 951 compounds (13.06%) with a predicted pIC₅₀ above 7.00, as seen in Figure 5(b). Applying criteria of high activity probability (probability ≥ 0.90) and a pIC₅₀≥7.00, we identified 68 potential FLT3 inhibitors for AML cells. Table S9 details these selected compounds. The scaffold analysis of the 68 compounds identified the most frequently occurring murcko scaffolds (Table S10), resulting in the selection of four compounds containing the most common murcko scaffolds for further molecular dynamics simulations and experimental validation, as shown in Figure 6.

Close

Figure 5.

(a) The chemical space distribution of 7280 molecules (b) The distribution of pIC₅₀ values for the 7280 molecules.

Export to PPT

Close

Figure 6.

The structural formulas of the four compounds selected for experimental validation.

Export to PPT

3.4. Molecular dynamics simulations

To evaluate the binding stability of the selected potential FLT3 inhibitors with FLT3 protein, we conducted a 100 ns molecular dynamics simulation on the docked complex of compound C1 and FLT3. Since the structural formulas of compounds C2-C4 are similar to C1, we will only discuss the case of compound C1 in this context. First, we examined the RMSD over the course of the 100 ns molecular dynamics simulation. RMSD is an essential measure that accurately represents the total displacement of all atoms from the reference conformation at any specific time point, acting as a primary indicator for evaluating the stability of the system. A detailed analysis of the data presented in Figure 7(a) reveals that as the simulation progresses, the RMSD values for both the complex structure and the protein exhibit a gradually stabilizing trend, with no significant fluctuations observed, indicating a progressively stable complex structure. Furthermore, to investigate the dynamic behavior of the small molecule on the protein surface, we focused on the interactions between the small molecule and the protein. By precisely calculating the distance between the centroid of the initial docking site residues and the centroid of the small molecule, as well as the distance between the ligand and the protein’s centroid, we enhanced our understanding of the ligand’s binding state to the protein. Figure 7(b) clearly illustrates that the distances between the ligand and the protein center, as well as between the ligand and its original binding site, remain highly stable without significant fluctuations. This finding strongly suggests that the small molecule consistently binds to the initial binding site on the protein throughout the simulation, demonstrating a stable binding relationship between the small molecule and the protein.

Close

Figure 7.

(a) RMSD and (b) Distance obtained from the MD simulation system between compound C1 and FLT3 protein.

Export to PPT

Considering the solvation effects, we integrated several parameters, including buried solvent-accessible surface area and interaction energy, to comprehensively assess the stability of the complex. Based on these considerations, we selected the stable-state trajectories of the complex and performed calculations using the Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) method. The binding energy and related energy terms obtained have been detailed in Table 3. In-depth analysis of the data in the table reveals that the van der Waals interaction energy (ΔE_vdw) within the complex is significantly higher than the electrostatic interaction energy (ΔE_ele), being specifically 2.6 times greater. Additionally, both are higher than the hydrophobic interaction energy (ΔE_nonpol). This indicates that van der Waals interactions dominate the binding energy composition, with electrostatic interactions playing a supplementary role, and hydrophobic interactions providing further support. Notably, the binding energy (ΔE_MMPBSA) between the small molecule and the protein reached -186.622±2.329 kJ/mol, a significant negative value indicating a very high binding energy and affinity between the small molecule and the protein. This finding provides important insights for further understanding and optimizing the interactions between the small molecule and the protein.

For a detailed structural and interaction analysis, we selected the conformation at the end of the simulation. As depicted in Figure 8, certain amino acid residues, including CYS828, ASP829, LYS644, GLU661, and CYS694, formed stable hydrogen bonds with the small molecule. Additionally, hydrophobic interactions such as π-Sulfur, π-π stacking, π-π T-shaped, and π-Alkyl were observed between the small molecule and amino acids like MET665, PHE691, PHE830, MET664, CYS828, LYS664, LEU818, ALA642, CYS694, and LEU616. Van der Waals forces were also noted between the small molecule and amino acids LEU668, GLU692, and GLY697. These results indicate that the small molecule binds securely to the protein, with high binding energy and affinity. Among the observed interactions, van der Waals forces are predominant, while electrostatic interactions are secondary, and hydrophobic interactions contribute significantly to the overall stability.

Close

Figure 8.

Binding interactions of compound C1 and FLT3: Interaction analysis of compound 1 with FLT3 based on the end of the MD simulations.

Export to PPT

3.5. Experimental validation: Biological activities assessment

After comprehensive screening using molecular docking, ML model predictions, and molecular dynamics simulations, we selected four compounds with potential inhibitory activity for assays on MV4-11 cell lines. In the experiments, we tested the effects of different concentrations of these compounds on the viability of MV4-11 cells and quantified their biological activity by calculating the IC₅₀. As shown in Figure 9, these compounds exhibited varying degrees of inhibition on the proliferation of MV4-11 cells at different concentrations. Among them, compounds C1-C4 all showed inhibitory effects on MV4-11 cells. Compound C2 demonstrated the strongest inhibitory effect with an IC₅₀ value of 1177 nM, indicating significant inhibitory activity against MV4-11 cells. Compound C2 followed with an IC₅₀ of 1731 nM. The other two compounds showed relatively weaker inhibition, with IC₅₀ values of 7280 nM and 8843 nM, respectively. The predicted values and experimental test values for these four compounds have been summarized in Table S11. These experimental results validate our screening method’s ability to effectively identify compounds with inhibitory effects on MV4-11 cell lines, further demonstrating the potential of integrating computational chemistry and ML-based screening strategies in early drug development. While the observed inhibitory activity is promising, the current IC₅₀ values have not yet reached the 100 nM threshold, indicating that further structural optimization through SAR (structure-activity relationship) exploration is needed. Key modifications could involve introducing halogen atoms (e.g., Cl, F) to enhance hydrophobic interactions with FLT3’s ATP-binding pocket, as well as optimizing hydrophobic substituents to strengthen van der Waals interactions.

Close

Figure 9.

Activity levels of compounds C1-4 in MV4-11 cells.

Export to PPT