Pharmacophore-based virtual screening, molecular docking, and molecular dynamics studies for the discovery of novel FLT3 inhibitors

Abstract FLT3 is considered a potential target of acute myeloid leukemia therapy. In this study, we applied a computer-aided methodology unifying molecular docking and pharmacophore screening to identify potent inhibitors against FLT3. To investigate the pharmacophore area and binding mechanism of FLT3, the reported co-crystallized Gilteritinib ligand was docked into the active site using Glide XP. Based on the docking results, we identified structure-based pharmacophore characteristics resistant to potent FLT3 inhibitors. The best hypothesis was corroborated using test and decoy sets, and the verified hypo was utilized to screen the chemical database. The hits from the pharmacophore-based screening were then screened again using a structure-based method that included molecular docking at various precisions; the selected molecules were further examined and refined using drug-like filters and ADMET analysis. Finally, two hits were picked out for molecular dynamic simulation. The results showed two hits were expected to have potent inhibitory activity and excellent ADMET characteristics, and they might be used as new leads in the development of FLT3 inhibitors. Communicated by Ramaswamy H. Sarma


Introduction
Acute myeloid leukemia (AML) is a heterogeneous category of leukemia that develops due to the clonal transformation of hematopoietic precursors with chromosomal rearrangements and numerous gene alterations (Saultz & Garzon, 2016). In AML, hematopoietic cells proliferate uncontrollably in the bone marrow. in an overabundance of aberrant white blood cells, infection, organ malfunction, and decreasing red blood cell production (Greim et al., 2014). According to GLOBOCAN, the world will have 474,519 new leukemia cases and 311,594 leukemia-related deaths in 2020 (https://gco.iarc.fr/). Approximately one-third of individuals with AML have a mutation in the FMS-like tyrosine kinase 3 (FLT3), caused by internal tandem duplications or a single mutation primarily affecting the tyrosine kinase domain (Daver et al., 2019).
FLT3 is a receptor tyrosine kinase (RTK) of class III found in normal and malignant lymphohematopoietic cells. It is crucial for the immunological response, hematopoietic stem cell formation, and proliferation (Meshinchi & Appelbaum, 2009). An external domain, a transmembrane helix, and an intracellular domain, which includes a juxtamembrane segment, a Tyrosine Kinase Domain (TKD), and a C terminal tail are all found in Class III RTKs. The FLT3 ligand binds to the transmembrane helix, causing dimerization and subsequent FLT3 autophosphorylation of tyrosine residues in the TKD, resulting in their activation and signal transduction through a variety of signaling pathways, including the RAS, SRC, and STAT5 pathways. AML has been shown to respond to several FLT3 inhibitors (Grafone et al., 2012).Several clinical trials have focused on different FLT3 inhibitors for each acute myeloid leukemia treatment setting -frontline therapy, treatment of relapsed refractory disease, and maintenance therapy after allogeneic HSCT (Sutamtewagul & Vigil, 2018). Clinical studies have been conducted on several inhibitors including Lestaurtinib, Midostaurin, Tandutinib, Crenolanib, Sorafenib, and Sunitinib (Wu et al., 2018). However, these inhibitors were not very selective for FLT3, and they also inhibited other class III RTKs such as KIT proto-oncogene receptor tyrosine kinase (KIT) and platelet-derived growth factor receptor (PDGFR) (Kiyoi, 2015). As a result, more powerful and selective FLT3 inhibitors such as Quizartinib, Crenolanib, Pacritinib, Ponatinib, and Gilteritinib were developed (Fathi & Levis, 2011;Tariq et al., 2021). So far, several AML drugs have been developed to target one or more faulty proteins. Last, on November 28, 2018, the FDA authorized Gilteritinib (XOSPATA, Astellas Pharma US Inc.) to treat adult patients with relapsed or refractory acute myeloid leukemia (AML) who have an FLT3 mutation (Smith, 2019;Stanchina et al., 2020). Despite the availability of these medicines, patients' overall survival rates remained at five years. To consolidate the progress made thus far in the hunt for drugs that can improve the overall survival of AML patients and overcome drug resistance of existing inhibitors (Reville & Kadia, 2020), the search for new inhibitors with fewer adverse effects is commendable and essential for life-saving treatments.
Present work combines ligand-based (molecular docking) and energy-optimized structure-based (e-pharmacophore) approaches for virtual screening of the free PubChem database. We used the PHASE module to create a structure-based pharmacophore model that offers the critical structural characteristics necessary for biological activity. The hits, as FLT3inhibitors, were recognized by utilizing molecular docking studies of pharmacophorematched compounds after the removal of pan-assay interference compounds (PAINS). Molecular dynamics simulations were carried out on selected poses to understand the dynamic behavior of the complexes and the stability of the protein-ligand binding Scheme 1. Our results will aid the further development of more powerful FLT3 inhibitors.

Software and computational analysis programs
In general, the various molecular modeling techniques used in computer-aided drug design play a very useful role in discovering new high-quality molecules more quickly and at lower cost than traditional methods that rely on manual screening in the laboratory.
In this computational study, we used Schrodinger's software to develop the predicted hypotheses based on e-Pharmacophore modeling procedures, structure-based virtual screening (SBVS), flexible ligand docking with the industryleading glide, 2 D ligand interaction diagrams, Prime MM/ GBSA, and advanced molecular dynamics simulations (MD) (Schr€ odinger Release, 2021a, 2021g). The use of this program is due to its high predictive power in discovering and improving new drug molecules faster and at a lower cost. In addition, the package of this software includes a number of useful applications for high resolution modeling, analysis, computation, and virtual screening. This is very encouraging for conducting in silico studies using the Schrodinger package without the need to look for other tools.
For visualization purposes, the BIOVIA Discovery Studio Visualizer 4.5.12 was used.

Pharmacophore hypothesis generation
Pharmacophores are three-dimensional representations of the characteristics of ligands essential to their biological function (Wilson & Lill, 2011). However, the methods applied here derive pharmacophore models by analyzing the interactions found in the binding sites of protein-ligand complexes, as opposed to ligand-based pharmacophore methods, which derive these arrangements solely from a series of active compounds. It is possible to create one or more pharmacophore hypotheses that describe the kind and spatial arrangement of the chemical characteristics necessary for a ligand to interact with the binding area's residues based on the ligandbinding site's characterization (Choudhury & Narahari Sastry, 2019). In order to have a more accurate and focused pharmacophore hypothesis, those elements crucial for ligand bioactivity should be chosen and included in the final model after being initially found by this method (Chandrasekaran et al., 2019).
The E-pharmacophore model of FLT3 was generated from the PDB ID 6JQR-Co-crystallized Gilteritinib ligand was docked in the generated grid using Glide XP docking. The protein-ligand coordinates the E-pharmacophore model was developed by using PHASE v.3 (Schr€ odinger Release, 2021b).
For this, the prepared protein-ligand complex was imported to the workspace, and default pharmacophore features such as hydrogen bond acceptor (A), hydrogen bond donor (D), aromatic ring (R), and hydrophobicity (H) were mapped. The best hypothesis was selected to screen compounds from PubChem database. The final hits were ranked based on the fitness score, ligand conformer that matches the hypothesis based on RMSD site matching, vector alignments and volume terms.
The E-pharmacophore model of FLT3 was generated using the crystal structure of FLT3 in a complex with Gilteritinib [PDB ID 6JQR]. -Co-crystallized Gilteritinib ligand was docked in the generated grid using Glide XP docking. The protein-ligand coordinates the E-pharmacophore model was developed by using PHASE v.3 (Schr€ odinger Release, 2021b).
For this, the prepared protein-ligand complex was imported to the workspace, and default pharmacophore features such as hydrogen bond acceptor (A), hydrogen bond donor (D), aromatic ring (R), and hydrophobicity (H) were mapped. The best hypothesis was selected to screen compounds from the PubChem database. The final hits were ranked based on the fitness score, ligand conformer that matches the hypothesis based on RMSD site matching, vector alignments, and volume terms.
The Generate DUDE Decoys program, located at http://dude. docking.org/generate, was used to create the decoy set. Researchers are able to create decoys for their assets using this tool, which is part of the Practical Decoy Database (Enhanced).
The decoy creation used the SMILES representation of the chosen inhibitor as input. The results of creating the decoy from the SMILES representation in the SDF 3 D structural format were converted using the Open Babel v2.3.2 program.
To identify inhibitors of FLT3, we have screened the PubChem database by searching compounds having more than 80% similarity instead of Gilteritinib compounds. The obtained compounds based on similarity sore were further filtered with Lipinski's Rule of 5 and reactive functional groups. The LigPrep protocols were used for ligands preparation by setting default parameters (Schr€ odinger Release, 2021c, Jorgensen & Tirado-Rives, 1988. The OPLS2005 force field was used for the minimization process (Kaminski et al., 2001). If provided in the input file, all chiral centers were kept; otherwise, up to 32 distinct stereoisomers per ligand were generated. For each compound, one low-energy ring conformation was produced. The ionization states and tautomer forms were measured using Epik at pH 7.0 ± 2.0.

E-pharmacophore based virtual screening
A pharmacophore-based virtual screening was carried out using the phase module of the Schrodinger suite to develop a subset of lead-like compounds with the requisite molecular characteristics for optimal binding to FLT3, as mapped by the E pharmacophore model Compounds had four sites on the e-pharmacophore hypothesis.
The final hits from screening were ordered according to their Phase fitness scores, which are based on vector alignments, volume terms, and RMSD site matching and assess how well the ligands match the chemical characteristics of the pharmacophore sites. The fitness rating is on a scale from 0 to 3, with higher ratings more comparable to Gilteritinib. The best hits were chosen based on fitness scores, and 820 compounds were selected for further studies.

Structure-based virtual screening (SBVS)
SBVS is a computational approach used in early-stage drug development campaigns to search a chemical compound library for novel bioactive compounds against a specific therapeutic target (Cheng et al., 2012). The initial screening utilized the e-pharmacophore model and identified 820 compounds, which shows the potential to interact favorably with the FLT3 active site. It the expected that the Pharmacophore model produces a lot of false positives in the search results; this is because they'll come across compounds with all of the characteristics but also a big side group that precludes the molecule from fitting into the active site. This structurebased screening was done on the selected ligands using Molecular Docking to determine this subgroup's most promising potential ligands. The GLIDE (Grid-based Ligand Docking with Energetics, version 7.8) module of Maestro 11.4 (Schr€ odinger Release, 2021e) was used for docking studies. Consequently, the ligands are sorted according to their affinity for the target, and the most promising molecules appear at the top of the list. Schema 1. Schematic representation of the methods followed in the current study

Molecular docking studies
The amino acid sequence of the protein FLT3 was reported from the UniProt database FLT3_HUMAN (code: P36888) at http://uniprot.org/ a freely available resource of protein sequences and functional information.
The X-ray crystallographic structures of the target protein FLT3 [PDB ID: 6JQR], obtained from the PDB database at https://www.rcsb.org in complex with Gilteritinib with a resolution of 2.20 Å were used for the study (Kawase et al., 2019).
The human FLT3 is a 993 amino acid long polypeptide with a molecular weight of 112.903 Da. It consists of three domains; (1) a 516-residue extracellular domain, (2) a 19amino acid transmembrane domain, and (3) a 429-residue cytoplasmic domain. The structure analysis was carried out using the PDBsum web server (http://www.ebi.ac.uk/thornton-srv/databases), which depicts the presence of 15 a-helices, six beta hairpins, four beta bulges, and 12 helix-helix ( Figure 1). The analysis revealed that Gilteritinib fits into the active conformation of FLT3 at the ATP-binding site (the 466 positions) in the activation loop. PDBsum online server was also used to check the validation of the FLT3 with the Ramachandran plot ( Figure 2), that's revealed that 93.5% residues in favorable regions, 6.5% residues in additional residue regions; 0.0% residues in generous regions; 0.9% residues in disallowed regions; Over all G-factor: 0.00.
The target protein was prepared using the protein preparation wizard in the Schr€ odinger suite (Schr€ odinger Release, 2021e). Bond orders, correction of formal charges to the hetero groups, and hydrogen atoms were added to all atoms in the system. At the same time, water molecules were removed, and the resulting structure was energy minimized until the heavy atoms converged to the RMSD of 0.30 Å using the OPLS_2005 force field (Shivakumar et al., 2012).
Glide's receptor grid generation module was used to create the energy grid for the Fms-Like Tyrosine Kinase 3 crystal structure Halgren et al., 2004). The docking grid box was specified by keeping the co-crystallized ligand Gilteritinib as the centroid. A grid was generated for Gilteritinib binding sites, i.e. (X: À 28.1226, Y: À 1.2759, Z: À 28.2055). To standardize the docking protocol, the co-  crystallized was docked to the corresponding active sites on FLT3 using a Glide docking module; the Extra precision (XP) docking technique is then used to identify the target-ligand interactions . This data was used for further development of the e-Pharmacophore model. Finally, the interaction of the co-crystal ligand with the active site residues of FLT3 protein was analyzed and cross-checked with Biovia Discovery Studio.

In-silico ADME prediction
After using the pharmacophore-based virtual screening and molecular docking studies, the top seven compounds were moved for ADME analysis. The Swiss ADME tool (Daina et al., 2017) was used for this. This tool assessed different pharmacokinetics parameters such as drug-like behavior, Lipinski and Veber's rules, bioavailability, synthetic accessibility, etc.
We have used the QikProp version 5.4 of the Maestro module (Schr€ odinger Release, 2021f) to predict the ADME properties of the best seven hits. This analysis includes aqueous solubility (Plog S), brain/blood partition coefficient (QP log BB), total solvent accessible surface area (SASA), log Khsa for human serum albumin binding (QPlogKhsa), octanol/water partition coefficient (QP log Po/w), predicted apparent MDCK cell permeability (QPMDCK), and human oral absorption

In silico toxicity studies
The toxicological evaluation of new drug entities is crucial in drug discovery and development. The ProTox platform (http://tox.charite.de/protox_II/) was used to evaluate the toxicity. This platform provides various substances' toxicity predictions, cytotoxicity, mutagenicity, carcinogenicity, immunotoxicity, and LD50 values.

Molecular dynamic and simulations
Molecular dynamics simulations were performed using the Desmond package for the free protein 6JQR and the complexes 6JQR-117675556 and 6JQR-155283200, which were selected based on the results of molecular docking simulations. All systems were solvated individually by putting them in an exhaustive water box of 10 Å with a single-point charge (SPC) (Jin et al., 2020;Wang et al., 2020). The force field OPLS3e was used for modeling the ligand-protein complexes (Choudhary et al., 2020;Gopinath & Kathiravan, 2021). The total charge of the systems subjected to the MD simulation was neutralized by adding sodium and chloride ions to the structure of the ligand-protein complexes. In addition, the energies of the systems were minimized to a minimum level with 2000 steps before running the MD simulation along a 100 ns trajectory in the NPT lattice (Al-Jumaili et al., 2021). In parallel, the ParticleMesh Ewald (PME) method was used to calculate the long-range electrostatic interactions while keeping a grid step of 0.8 Å (Rajagopal et al., 2021). The Nose-Hoover thermal algorithm and the Martina-Tobias-Klein method were used to generate slow heating of the systems with fewer than 300 K and 1 bar pressure (Martyna et al., 1992). The Desmond package simulation interaction diagram tool analyzes the thermodynamic parameters of the systems investigated in the aquatic environment (Daoui et al., 2022a(Daoui et al., , 2022b(Daoui et al., , 2022c. This is done by generating the timing of variations in each system's parameters (Total and potential free energies, temperature, pressure, and volume). Furthermore, the energy parameters generated by the MM-GBSA (Molecular mechanics-generalized born surface area) simulation were predicted by the mmgbsa_thermal.py algorithm provided in the Schrodinger 2020-3 software package. This was done to predict the stabilization energy levels of potential interactions between the selected ligands and the 6JQR target receptor.

Pharmacophore generation
In the present study, an e-pharmacophore hypothesis was developed to screen the FLT3 protein inhibitors by using its crystal structure in complex with a potent, broad-range noncovalent inhibitor (the native ligand) after XP docking. The details of this step are explained in the section (Materials and methods). The generated e-pharmacophore model contains four-featured, i.e. AHRR [one hydrogen bond acceptor (A), one hydrophobic (H), and two aromatic rings (R)]. The pharmacophore concept of AHRR in the planar form is provided in Figure 3. The hydrophobic site H14 lies on the C14 atom, whereas the acceptor, A5, lies on C ¼ O (O5 atom), and groups R18 and R19 lie on the rings of the ligand.

E-pharmacophore validation
The generated pharmacophore model was validated to see if it can accurately predict the activity of novel compounds discovered through database screening or created from scratch.  Validation of a pharmacophore model is a critical step before it can actually be used for virtual screening (Kaserer et al., 2015). The model was employed to screen a testing set database that included 21 known FLT3 inhibitors with experimental activity and 50 inactive molecules generated from the DUDE Decoys tool. For validation, various statistical parameters, including EF, RIE, BEDROC, AUAC, and ROC, were calculated that are presented in Table S1-S4. The values of enrichment factor (EF) to recognize active compounds. Validation of the FLT3 e-pharmacophore hypothesis revealed that EF in the top 1% of the decoy dataset is 27%, demonstrating that FLT3-specific e-pharmacophore is 27-fold efficient in detecting true positives/actives from the entire dataset. In the case of ROC, values closer to 1 correspond to a reliable pharmacophore in differentiating active compounds from non-active compounds, whereas RIE higher than10.0 are generally preferred for that purpose. ROC score, RIE, and AUAC values were calculated as 1, 14.22, and 0.98, respectively. Thus, FLT3 e-pharmacophore was statistically significant in picking the actives from the decoy dataset. Statistical significance of FLT3 -specific e-pharmacophore was also validated by calculating BEDROC. Contrary to EF, which has been detailed above, BEDROC seeks to measure the early enrichment of the actives. BEDROC values were calculated at different tuning parameter values (a ¼ 8.0, a ¼ 20.0, and a ¼ 160.9) and found to be 0.999, 0.999, and 1, respectively.
The receiver operating characteristic (ROC) curve of the pharmacophore model is shown in Figure S1. Moreover, the ROC value of 1 indicates that the pharmacophore model has a probability of 1 to effectively screen active molecules with required inhibitory activity (Hajian-Tilaki, 2013). It is also apparent from Figure S1 that the ROC curve initially progresses vertically, implying the model's effectiveness in finding and ranking the actives at the beginning of the screening process.

Structure-based docking
Initially, the pharmacophore model was applied to search the ligand database. The hit compounds were considered for a structure-based molecular docking study to observe the binding interactions between ligand and catalytic residues at the active site of the selective receptor to choose the best hits with the best glide docking score. The 820 compounds were still large to put for docking in extra precision mode (XP), so docking in standard precision (SP) mode was carried out to refine the selection to a maximum of 40 compounds. SP docking requires less CPU (central processing unit) time, and its scoring function includes lesser terms than the XP scoring function, which uses explicit-water technology and descriptors. The XP docking approach effectively eliminates false positives and improves the connection between excellent postures and good results. The sampling approach of XP mode is based on an anchor and improved growth strategy, and it is less forgiving than the SP mode.

Interaction profiles of hit molecules with FLT3 RECEPTOR
Finally, in the last stage of the VS, we selected top hits compounds on the basis of XP docking score of Gilteritinib < À 8.8 kcal/mol at the active site residues of receptor tyrosine kinase (We found only seven compounds, CID 155283200, CID152889583, CID145613128, CID117675549, CID117675556, CID117675577, and CID117675745) having the extra precision glide docking scores less than that of the native ligand hits ranging from À 8.065 to À 9.764 as provided in Table 1. The hydrogen bond interactions, Van der Walls interactions, electrostatic interactions, aromatic interactions, and Ligand-binding interactions were also studied According to the result of molecular docking, the reference ligand (Gilteritinib) was determined to have a binding energy of À 8,038 kcal/mol. It can be inferred from the interaction diagram that the oxygen atom of Gilteritinib forms a conventional hydrogen bond with CYS694 residue, and the amino group of Gilteritinib also was found to form two conventional hydrogen bonds with Glu692, and Phe830. We also see the hydrophobic interaction with residues Leu616, Leu818 type pi_alkyl lie at the ring pyrazine ( Figure S2). Recent literature also shows that these amino acid residues are present in the binding pocket of the FLT3 receptor   the binding poses disclosed that the compounds oriented in the FLT3 binding cavity similarly bind with residues from both the membrane-binding and catalytic domains. The mapping of the compounds onto the pharmacophore and the orientation of FLT3 inhibitors in the active site with interacting residues in Figure S3A. Figure S3B shows that the seven hits have the complete mapping with the four features of pharmacophore model AHRR and the better alignment of compound CID 155283200 with features, resulting in a higher fit score (Table 1) and forming stronger interactions with ATP of FLT3 receptor. In comparing the binding interactions of proposed molecules with the co-crystal ligands, it was found that all catalytic amino residues for co-crystal ligands actively interacted with the proposed molecules; the interaction patterns of the hit molecules were also examined. Interestingly, the hit molecules CID 155283200, CID152889583, CID145613128, CID117675549, CID117675556, CID117675577, and CID117675745 were also forming hydrogen bonding interactions with Glu692 and Cys692 hydrophobic interaction with Leu616, Leu818, Val675,   Figure S3 by 2 D interaction diagrams. The best performing protein-ligand single-docked coordinate of each complex was further taken for the MD simulation study.

Compliance with a standard range of drug-likeness, bioavailability, synthetic accessibility, and alerts for PAINS
The proportion of global pharmacokinetics-related failure in the later stages of drug development can be reduced by applying the in-silico drug-likeness study during the early phases of drug discovery. This study helps to minimize the drug development cost (Schneider, 2013). In this study, we estimated the Drug-likeness profiling of hit compounds with SwissADME (http://www.swissadme.ch) and provided the result in Table 2. The physicochemical parameters such as PM, TPSA, HBA, HBD, LogP, and the number of rotative bonds were calculated for hit molecules and compared with known drugs. All the compounds were also screened with Lipinski's rule of five, resulting in Table 2. The result shows that all seven hit compounds met Lipinski's rule, indicating they had drug-like potential. The screening process with Veber rules demonstrated that all compounds meet the criteria of drug-likeness assessment except compound CID152889583, which had one violation (TPSA >140). In drug development, the bioavailability of an orally administered drug, or the percentage of the dosage that reaches the circulation, is essential. The Abbot Bioavailability score predicts the chance of the compound to have at least 10% oral bioavailability in rats. The candidate compounds showed a score of 55%, indicating good bioavailability. Along with the ADMET study, the SwissADME program was also used for Pan-assay interference compounds (PAINS). The PAINS method is used to identify potentially problematic fragments that yield false-positive biological output (Baell & Holloway, 2010). The results of this study indicate that compounds 155283200 and 117675556 did not show any such fragment. The rest of the compounds show violations because the query dataset's chemical structure contains fragments such as aniline. The compounds were also analyzed for synthetic accessibility, and the results revealed that the compounds could be synthesized easily (Synthetic accessibility <3). The BOILED-Egg graph (Daina & Zoete, 2016) for compounds analyzed in this study is shown in Figure S4. The charts plotted between WLOGP and TPSA (topological polar surface area) values. This graph shows that only one molecule is out of range, i.e. CID145613128, and other molecules have shown GI absorption value by occupying the white space on the BOILED-Egg chart.

Prediction of ADMET properties
The QikProp module of Schr€ odinger was used to predict the ADME properties of the seven hit molecules. The results of the ADME properties are shown in Table 3. Jorgensen's rule of three was applied to determine the bioavailability of each compound by calculating its solubility, permeability, and liver first-pass metabolism through the following rules:-logS (> À 5.7), PCaco (>22 nm/s), and primary metabolites (PM) (<7). These rules violations are essential for optimizing biologically active compounds and should not be more than 1. Table 3 lists Jorgensen's rule result, and all compounds showed ADME properties in an acceptable range. The #star descriptor was also considered, which indicates how many characteristics of each chemical deviate from the recommended range. As a result, compounds with a lower #star value have more incredible drug-like qualities. The seven hit molecules were found to have a #star value of 0 as that of the known drug.
Calculations related to serum protein binding, blood-brain barrier (QlogBB and apparent MDCK cell permeability), gut blood barrier (QPPCaco), predicted central nervous system activity (CNS), skin permeability (QPLogKP), and human oral absorption indicated that these values for the hit compounds with the standard ranges generally observed for the drug.
Toxicity is one of the last characteristics of the molecule's ADME/Tox study. We employed the web-servers ProTox-II to predict the toxicity parameters. The results suggested that none of the compounds shows hepatotoxicity, carcinogenicity, immunotoxicity, mutagenicity, and cytotoxicity except compound CID 152889583 and CID 117675549, which shows immunotoxicity, and compound CID 145613128, which shows carcinogenicity.
Regarding acute oral toxicity, based on the ProTox results, CID 155283200, CID 145613128, CID 117675549, CID 117675577, and CID 117675745 had LD50 values between (300 < LD50 � 2000) and in-class IV Global Harmoni System (GHS) indicating that it could be harmful if swallowed. The two other compounds, CID 152889583 and CID 117675556 have LD50 values between (2000 < LD50 � 5000) and in-class V GHS, so they may be harmful if swallowed (Table 4).

Molecular dynamics study
The dynamics of the 6JQR protein in its forms, free and complexed with ligands CID 117675556 and CID 155283200 under an aqueous environment were investigated. This was done to evaluate the stability level of the studied systems as a function of MD simulation time of 100 ns. For this purpose, we used the complexes 6JQR-CID 117675556 and 6JQR-CID   155283200 obtained through molecular docking studies. Figure 4(A-C) shows the time frame of RMSD and RMSF index variations for the examined systems.

RMSD (root mean square deviation) and RMSF (root mean square fluctuation) analysis
Molecular dynamics simulations were performed over a 100 ns time trajectory to understand the underlying molecular insights into the binding of ligands 117675556 and 155283200 in the 6JQR protein pocket. The MD simulation trajectories were initially analyzed by calculating the RMSD of each system. The RMSD values of the free 6JQR backbone were adopted as a reference to fit the RMSD values of 6JQR-117675556 and 6JQR-155283200 complexes. From Figure 4A, we can see that the backbone atoms of 6JQR protein in its free and complex forms remain well stable throughout the 100 ns simulation time frame. Some slight deviations were observed in RMSD values of amino acid residue atoms of 6JQR free protein. So that RMSD values were observed to increase from 1.64 Å at 13 ns to 2.38 Å at 50 ns before stabilizing at 55 ns with a value of 2.07 Å until the end of MD trajectory. In addition, we can observe that RMSD values of backbone atoms in complexes 6JQR-117675556 and 6JQR-155283200 remained stable throughout the MD simulation. The presence of some small deviations of RMSD values (less than 2.28 Å) was observed in complex 6JQR-155283200 at 50 ns.
The average RMSD values for the free 6JQR backbone complexed with ligands 117675556 and 155283200 did not exceed 1.94 Å, 168 Å, 1.62 Å during the 100 ns simulation run. The RMSD values observed for all three systems did not reach the 3 Å thresholds, which explains the high potential stability of ligands 117675556 and 155283200 in the active pocket of the 6JQR protein. RMSF computations for the backbone atoms of 6JQR free of and complexed with the ligands (117675556, 155283200) were analyzed ( Figure 4B) The RMSF values obtained indicate that all residues of free and complexed 6JQR with ligands 117675556 and 155283200 are less than 2.68 Å. This confirms the high stability of the 6JQR protein residues in their free conformation and the conformation complexed with the selected ligands 117675556 and 155283200. Also, the RMSF values of the atoms of ligands 117675556 and 155283200 were computed to describe the dynamic response of the ligands within the 6JQR protein pocket. The data presented in Figure 4C indicate that some fluctuations in the ligand structures did not exceed the thresholds of 1.55 Å and 1.90 Å for the two ligands, i.e. 117675556 and 155283200. The fluctuations observed in the ligand structures may explain that the structures of molecules 117675556 and 155283200 have undergone numerous dynamic mutations at their binding sites in the 6JQR protein pocket. The observed fluctuations in the structures of ligands can be explained by the flexible structure character of the examined ligands. Despite these fluctuations, the molecular structures of the two ligands remain stable, with the 6JQR backbone atoms over 100 ns in aqueous conditions.

Ligand-protein contacts dynamics
Figures 4 and 5 show the timescale of the contact dynamics between 6JQR protein and ligands (117675556 and 155283200).
Visualizing the protein-ligand complex will help to provide better insight into the MD simulation results. From Figure  5(A-C), we can observe that the active residues Glu692, Cys694, and Phe830 are among the most important residues observed to play an essential role in maintaining the stability of ligand 117675556 in the 6JQR protein pocket. Glu692 provided hydrogen bonding to the nitrogen atom for approximately 100% simulation time. Similarly, Cys694 contributed hydrogen bonding to the oxygen atom of the group amide for more than 100% of the time of the MD simulation, indicating the high stability of the ligand 117675556 bonds with the active site of the Cys694 residue.
Similarly, the side chain of residue Phe830 contributed another hydrogen bond to the nitrogen atom up to 70% of the simulation time of 100 ns. When the contact strength decreased to less than 40%, temporary bridging interactions were observed between water and residues Cys695, Tyr696, Asn816, and Asp698. Some hydrophobic interactions were also observed between the ligand 117675556 structure and residues Val624, Tyr693, Leu818, Cys694, and Phe830, with percentages ranging from 10% to 65% of the MD simulation path.
From Figure 6(A-C), we can notice that the active residues Leu 616, Glu692, and Cys694 are the most important residues whose substantial contribution is maintaining the ligand's stability 155283200 in the 6JQR pocket. Leu 616 and Glu692 gave two hydrogen bonds to the nitrogen and oxygen atoms of the azanide group at 75% and 100% of the MD pathway, respectively. When the contact strength decreased to less than 20%, temporary bridging interactions were observed between water and the Asp698 residue. Some hydrophobic interactions were also observed between the ligand CID 155283200 structure and residues Leu 818, Tyr693, Phe830, Val624, and Ala642, with average rates between 15% and 65% of the MD simulation pathway time.

Properties analysis of ligands 117675556 and 155283200
To further validate the high stability hypotheses of the examined complexes in the aqueous environment, we predict the most important dynamic properties of ligands 117675556 and 155283200, contributing to the high interaction stability of ligands 117675556 and 155283200 to 6JQR protein active sites. Figure 7(A,B) shows the summary of the six examined properties (Ligand RMSD, Radius of Gyration (rGyr), Molecular Surface Area (MolSA), Intramolecular Hydrogen Bonds (intraHB), Solvent Accessible Surface Area (SASA), Polar Surface Area (PSA)) calculated for the two ligands in the 6JQR-117675556 and 6JQR À 155283200 complexes, respectively. From Figure 7(A,B), we can comprehend that the RMSD values of ligand 117675556 remained generally stable between 0.4 Å and 1.2 Å throughout the MD simulation (run set at 100 ns). From Figure 5B, we can also realize that the RMSD values of ligand 155283200 remained stable between 0.5 Å and 1.6 Å at 76 ns, after which the RMSD values increased slightly to 1.8 Å at 80 ns. After that, the RMSD values decrease again to stabilize at 0.6 at the end of the simulation. Generally, the slight variations in the RMSD values for the ligands (<3 Å) during the MD simulation run (100 ns) explain the high stability of the ligands 117675556 and 155283200 scaffolds in the 6JQR protein pocket. The average stability of rGyr values in the 100 ns path for the MD simulation time for ligands 117675556 and 155283200 is 4.80 Å and 4.39 Å, respectively. With some slight fluctuations in the radius values of ligand 117675556 between 4.59 Å and 4.84 Å at 16 ns to 30 ns and 70 ns to 81 ns. After 81 ns, the stability of the rGyr values of ligand 117675556 in the 6JQR pocket was observed to be approximately 4.6 Å until the 100 ns duration of MD simulation. For ligand 155283200, we can also notice that there are slight fluctuations not exceeding 4.8 at 76 ns to 80 ns, then the rGyr values of ligand 155283200 remain stable at 4.35 Å until 100 ns. The intraHB, MolSA, SASA, and PSA plots presented in Figure 7 The analyses of complex 6JQR-117675556, complex 6JQR-155283200, and free protein 6JQR system dynamics indicate the high stability of the ligands in the 6JQR protein pocket. To confirm this hypothesis, the evolution of physical and chemical parameters thermodynamics such as the total and potential energies (ET, EP), temperature (T), pressure, and volume (V) of 6JQR complexes were evaluated and provided in Figure 8. From Figure 8A, we can notice that the parameters ET, EP, T, P, and V of the 6JQR-117675556 complex remained constant during the whole simulation period at these average values of À 87543.589 KJ/mol, À 107198.608 KJ/mol, 298.698 K, 1.316 bar, 319680,51 Å3, respectively. Figure 8B also shows that the parameters ET, EP, T, P, and V of 6JQR-155283200 remained constant throughout the simulation period, with mean values of À 85848.378 KJ/mol, À 106063 KJ/mol, 298.706 K, 1.332 bar, and 329263.042 Å3, respectively. The stability of the parameters ET, EP, T, P, and V throughout the simulation period validated the stability status of the investigated 6JQR-117675556 and 6JQR-155283200 complexes.

Evaluation of MM-GBSA parameters
We performed MM-GBSA simulations to evaluate the parameters of free binding energies between the ligands (CID 117675556 and CID 155283200) and target protein receptor tyrosine kinase (RTK) of class III (6JQR). This procedure is the final step in checking the stability levels of the examined systems in the aquatic environment. Table 5 shows the average values of calculated MM-GBSA energy parameters (total binding free energies of Coulomb energy (Coulomb), covalent bonding (Covalent), hydrogen bonding (H-bond), lipophilic bonding (Lipo), p-p stacking interaction (Packing), solvent generalized bonding (SolvGB), and Van der Waals bonding energy (VDW)). Figure S6 illustrates the average values of the calculated energy parameters of the examined systems over the 100 ns run of the MD simulation. From Table 5, we can observe a strong convergence between the parameter scores of the binding energies relative to the 6JQR-117675556 and 6JQR-155283200 complexes, which can be confirmed by the coinciding distribution of these parameters in the two complexes, i.e. 6JQR-117675556 and 6JQR-155283200 ( Figure S6). Overall dynamic parameters obtained via MD simulations indicate the high stability of CID 117675556 and CID 155283200 ligands in the pocket of the target protein of 6JQR. Thus, the inhibitory potential of CID 117675556 and CID 155283200 against FLT3 is increased. Therefore, the structures of these two compounds should be suitable for treating acute myeloid leukemia by targeting FLT3.

Conclusion
The present study employed a structure-based pharmacophore and structure-based virtual screening to identify potential FLT3 inhibitors from the PubChem library. The developed e-pharmacophore model consisted of one hydrogen bond, one hydrophobicity, and two aromatic groups. Extensive enrichment analysis was used to assess the model's quality.
Seven compounds with the diverse structure were obtained based on high docking scores and favorable docking patterns. These hits were further analyzed for their pharmacokinetic properties using ADME analysis. After a series of screening, two compounds were selected for molecular dynamics simulation studies to determine the stability of their binding in the active site of the FLT3 receptor, and the results showed that the CID 155283200 and CID 117675556 bind strongly to the enzyme active site of FLT3 receptor and hence they could be purposed acute myeloid leukemia against.
Overall, this research shows that combining virtual screening molecular dynamics simulations can help researchers investigate AML-targeted small molecule screening.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
The author(s) reported there is no funding associated with the work featured in this article.