Intermediate-based virtual screening of c-Kit kinase inhibitors as potential anti-tumor agents via ab inito folding, molecular dynamics simulation, and molecular docking

2.1 Secondary structure profiling

The secondary structure profile of the activation loop was generated with PEP-FOLD 4.0 application(Rey et al., 2023) via RPBS Web Portal (https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::PEP-FOLD4). The sequence of the activation loop of KIT kinase was used for secondary structure profile prediction, and the detailed sequence used for this experiment was “THGRITKICDFGLARDIKNDSNYVVKG”.

2.2 Ab inito folding of kinase domain within c-Kit kinase

The crystal structure describing the autoinhibitory conformation (PDBID:1T46) and active comconformation (PDBID:1PKG) of the tyrosine kinase domain within c-Kit was retrieved from Protein Data Bank. After removing water and co-crystallized ligand, the structural model was imported into the PyRosetta. The folding protocol was constructed according to the literature with adjustments. Briefly, a loop region between Cys 809 to Asn 831 was defined, and the centroid conformation of the defined loop region was generated with FragmentMover and CCDLoopClosureMover. After recovering the side chain conformation with PackRotamersMover, the full atom structure was subjected to the final refinement with LoopMover_Refine_CCD. The “score3” and “talaris2013” energy score function from PyRosetta was employed to evaluate the potential energy of the centroid and full atom structural model, respectively. For each initial structure, 5000 conformations (decoys) were generated. The root-mean-square deviations (RMSD) of the A-loop compared to the two reference structures were calculated with PyMOL, and the K-means clustering was performed by the scikit-learn package(Pedregosa et al., 2011) with Python. The three-dimensional structure of the generated decoys was visualized by PyMOL 2.5 (Open-source version).

2.3 Molecular dynamics simulation

The molecular dynamics (MD) simulations were performed to obtain the stable intermediate conformation. The hydrogen atoms were removed, and the protein structures were parameterized with the Amber 14SB force field(Maier et al., 2015). Then, the protein was embedded into a cubic box and dissolved within the TIP3P water model. The total atom number of the prepared systems was around 35,600, and the initial dimension was around 70.6  Å × 78.7  Å × 73.5  Å.

To remove the defects in structural models and to resolve the collision between hydrogen atoms, a three-staged energy minimization protocol was employed as a previous study. Briefly, the complex was subjected to an all-heavy atom position-restrained energy minimization (lasted for 500 steps), backbone atoms position-restrained energy minimization (for 500 steps), and an energy-minimized without restraints for another 1,000 steps. The energy-minimized system was gradually heated to 310  K from 0 K in 100 ps with all heavy atoms restrained. Finally, the restraints were gradually removed in a constant pressure environment (NPT) for a total of 2.5  ns. The simulation time step was set to 2 fs, and the van der Waals interactions were calculated at every step. The long-range electrostatic interactions were calculated every 5 steps by the particle mesh Ewald (PME) method(Fukuda et al., 2014) with a grid size of about 1  Å and a tolerance of 10⁻⁶. The cut-off for the van der Waals interactions was set to 9  Å. The temperature and pressure were controlled using Langevin dynamics and the Langevin piston barostat method. The SHAKE method was used to restrain hydrogen atoms and the tolerance was set to 10⁻⁸. Atomic coordinates of all atoms were saved every 1 ps. The production NPT run was performed at 310  K and 1 atm for at least 100  ns for each intermediate conformation by NAMD 2(Phillips et al., 2020). After obtaining the trajectory of the NPT production run, the RMSD, root-mean-square-fluctuations (RMSF) corresponding to the reference structure, as well as the secondary structure profile were calculated with Gromacs. The secondary structure of the protein was calculated by DSSP (Joosten et al., 2011) and visualized with do_dssp in Gromacs.

2.4 Preparation of compounds dataset and database

To evaluate the descriptive power of the structural model, a small dataset containing 6 known KIT inhibitors, including imatinib, masitinib, bezuclastinib, avapritinib, BLU-263, and elenestinib was constructed according to literature(Akin et al., 2022). To identify the potential inhibitors against KIT with D816V mutation, compounds from the public database of available kinase inhibitors currently in clinical trials (from phase I to IV), PKIDB(Carles et al., 2018, Bournez et al., 2020), were retrieved for virtual screening. Then, the three-dimensional structure of compounds in the test dataset and PKIDB were generated with the Ligprep application in Schroginder Suit. The bond angles and bond ordering were assigned, and the minimization using the OPLS-2005 force field with default settings.

2.5 Molecular docking

The protein structure of c-Kit obtained from molecular dynamics simulation was processed with the Protein Preparation Wizard application of Schrödinger Suite (Schrödinger 2020). The ionization state of the protein was calculated by PROPKA application(Olsson et al., 2011), and the hydrogen atoms were added to optimize the hydrogen bonds. Finally, the optimized structure was parameterized with the OPLS-2005 (Optimized Potentials for Liquid Simulations) force field and the energy minimization was performed with default settings(Sastry et al., 2013). The Glide application of the Schrodinger suit was employed for molecular docking. The grid which centered on the DFG motif within the activation loop was generated with a default value was generated using the default parameter, and the binding free energy was predicted by SP (standard precision) mode in the Glide application.

2.6 In vitro cytotoxicity test

The HMC1.1 cells containing V560G in the juxtamembrane domain of KIT and imatinib-resistant HMC-1.2 cells containing V560G and D816V dual mutations were maintained in Iscove's modified Dulbecco's medium (Gibico, Guangzhou, China) supplemented with 10 % fetal calf serum (Gibico, Guangzhou, China) at 37 °C, 5 % CO₂. After treatment with Vorolanib of the indicated concentration for 24 h, the cell viability was measured with the Cell Counting Kit − 8 assay according to the literature(Yang et al., 2021).

2.7 Statistical analysis

Pooled data were analyzed by GraphPad prism software and visualized as the mean ± SEM. Differences between each group were analyzed using one-way ANOVA followed by Tukey post-hoc test. Probabilities (p) less than 0.05 were considered to be significant.

3.1 Secondary structure prediction revealed the presence of short α-helices adjacent to the DFG motif was one of the characteristics of the intermediate state

By comparing the crystal structure describing the kinase domain of c-Kit in autoinhibitory and activated states, we found a remarkable conformational switch of the A-loop. In its autoinhibited state, the A-loop adopted a “closed and folded” state, which characterized the presence of a short helix between C810 and G812, as well as a beta strand between D820 and L831(Mol et al., 2004). In the active state, the activation loop and DFG motif adopted an extended conformation, and the helix and strands were replaced by random coil(Mol et al., 2003) (As shown in Fig. 1A and Fig. 1B). As accumulating evidence demonstrated the activation of c-Kit kinase, as well as the member of SRC kinases required secondary structure changes of the activation loop(Ozkirimli and Post 2006), we hypothesize the metastable state of the c-Kit kinase domain might adopt a different secondary structure profile from identified crystal structures.

Close

Fig. 1

Secondary structure prediction revealed the presence of short α-helices adjacent to the DFG motif was one of the characteristics of the intermediate state. The crystal structure describing the KIT kinase in auto-inhibitory and activated conformation was shown in A and B, respectively. C. The results from the secondary structure prediction of residues in the activation loop, the presence of a-helices, b-strands, and the coil were colored red, blue, and green, respectively, while the representative structure was illustrated in D.

Export to PPT

To obtain more information about the secondary structure profile about the possible metastable state, we performed a secondary structure prediction towards the activation loop. The results identified the presence of a β-sheet-containing structure within the D820 to L831, as well as a short β-strand located close to the beginning of the DFG motif within the kinase domain (Fig. 1C), which was in line with the available crystal structures. Interestingly, the secondary structure prediction revealed the residues within F811 to I817 adopted α-helices (colored in red) and random coil (colored in blue) with similar possibilities. However, the presence of short α-helices ranging from F811 to I817 (shown in Fig. 1D) was not observed in all available crystal structures. One the other hand, the presence of short α-helices was observed in the metastases found in many protein kinases, such as the anaplastic lymphoma kinase(Patil et al., 2021), Bruton tyrosine kinase(Sultan et al., 2017), Src kinase(Shukla et al., 2014) and Abl kinase(Paul et al., 2020). As a result, we consider the short helix within A-loop might be also a characteristics of metastable state in c-Kit enzyme, and we used it as a criterion for the identification of possible metastable states.

3.2 Ab initio folding and K-means clustering identified nine potential metastable conformations of c-Kit kinase, which showed remarkably different secondary structure profiles

Several computation methods were employed to study the dynamic properties of the kinase domain within c-Kit and identify the possible metastable state of c-Kit, and the results demonstrated the conformational alterations of residues within JM domain, P-loop, as well as the random coil located in N-lobe was closely related to the conformational switch between inactive and active state(Yang et al., 2009, Laine et al., 2011, Shukla et al., 2014, Paul et al., 2020). However, the detailed conformational information about the A-loop, which played an essential role in the activation process, was still limited. One of the possible reasons was the insufficient conformational sampling made by molecular dynamics simulation and enhanced sampling method(Del Alamo et al., 2022, Zhu et al., 2023). As a result, the PyRosetta package, which incorporated the Monte Carlo sampling technique was employed. Rosetta and PyRosetta were designed for de novo protein structure prediction, assembling fragments, which provided a promising platform to accomplish an extensive and sufficient sampling in the protein conformational space(Heilmann et al., 2020). Accordingly, Rosetta and PyRosetta were already widely used for the modeling and design of anti-bodies(Schoeder et al., 2021), enzymes(Bhattacharya and Cheng 2015), and other macromolecules(Leman et al., 2020).

To expand the conformational space of the activation loop within the kinase domain, and identify the structural mode of the intermediate, an ab inito folding-based structure predicting process was employed. As illustrated in Fig. 2A, the crystal structure describing the autoinhibitory and active state were employed as starting conformation, and subjected to ab initio folding with PyRosetta. After loop closure, side-chain packing, and energy minimization, the predicted conformation (decoys) was generated. To maximize the conformation pool, 5000 decoys for each starting structure were generated. Then, the RMSD of each decoy compared to the starting structure was calculated and plotted in Fig. 2B. Based on the result from the Silhouette analysis, K-means clustering of RMSD scatter identified nine clusters based on the RMSD plot, and the decoys belonging to each cluster were ranked based on their predicted energy. The top-scored poses bearing the lowest potential energy were selected as representatives, and visualized in Fig. 1C. Based on the conformation of the DFG motif, the identified metastable conformation could be subdivided into two different groups. In the first group (conformation a to d), the A-loop adopted a closed and folded structure. The β-sheet ranging from N822 to L831 was preserved, and the A-loop resided in the cleft between αC helices. Similar to the KIT kinase in its auto-inhibitory form, the DFG motif preserved its DFG-out conformation. While in conformation e to i, the DFG motif adopted a DFG-in conformation, which was similar to active KIT kinase. Besides, we noticed a remarkably different secondary structure profile within the predicted metastable intermediates. Based on the results from secondary structure prediction, we found that conformation b and conformation e, which possess the short α-helix within the A-loop were the structural model of potential metastable conformation.

Close

Fig. 2

Ab initio folding and K-means clustering identified nine potential metastable conformations of c-Kit kinase, which showed remarkably different secondary structure profiles. A. The illustration of the intermediate generation process. The crystal structure of the kinase domain of c-Kit in the autoinhibitory state (PDBID: 1 T46, colored cyan) and active state (PDBID: 1PKG, colored pink) were used as the initial structure for ab inito folding. For each initial structure, 5000 decoys were generated. The RMSD of each decoy related to the initial structures was plotted in B, and the representative structure of each cluster was visualized with a cartoon in C..

Export to PPT

A careful analysis of the identified metastable state revealed the unfolding of the A-loop happened before DFG-flip, indicating the unfolding of the A-loop was one of the pre-requests for DFG-flip and kinase activation. Besides the identified metastable states b and e which may possess higher stability, we also found four intermediate conformation f to the conformation i, which was similar to the “DFG-up”, “DFG-in” and “DFG-down” states identified in a metadynamics simulation(Jiang et al., 2022).

3.3 Molecular dynamic simulation identified one possible metastable state with promising thermal stability

In vitro and in silico studies demonstrated there were a variety of intermediates that existed between the active and inactive states, however, their thermal stability differed significantly(Laine et al., 2011). Accumulating evidence revealed that only the stable intermediate along the transition pathway may play an essential role in regulating the catalytic activity or molecular function(Kannan and Zacharias 2014, Ghosh and Ranjan 2020). The formation of stable intermediate was accompanied by the exposure of buried druggable pockets, which made them an attractive target in the discovery of novel kinase inhibitors (Umezawa and Kii 2021). Accordingly, the thermal stability of two potential intermediates was investigated with molecular dynamics simulations, and the secondary structure profile, as well as the RMSD of the A-loop referring to the initial structure, was measured. The higher stability of the secondary structure profile and the lower RMSD value indicated the increased thermal stability of the simulation system.

As shown in Fig. 3, the molecular dynamics simulation did not alter the overall secondary structure of KIT kinase, however, the α-helical structure of A-loop within the conformation b was gradually replaced with 3-helix (Fig. 3A, middle panel), indicating the helix was less stable. While, the α-helix within the conformation e was preserved during the simulation, which indicated the higher stability compared with conformation b (Fig. 3B, middle panel). On the other hand, we observed less variation in the RMSD plot at 80 to 100 ns within conformation e, further indicating the increased thermal stability of α-helix within the activation loop (Fig. 3A and B, lower panel).

Close

Fig. 3

Molecular dynamic simulation identified one possible metastable state with promising thermal stability. The secondary structure profile of conformation b and conformation e was obtained from molecular dynamics simulation and was shown in the upper panels of A and B. The secondary structure of the A-loop (residues between 180 and 190) was plotted in the middle panel of A and B. The RMSD of the A-loop along with the simulation was shown in the lower panel of A and B.

Export to PPT

In addition, the RMSF of two structural models was calculated. Comparing the RMSF profile of two structural model, we observed the overall RMSF, as well as the RMSF of A-loop (high-lighted in red frame), was smaller within confirmation e than conformation e (As shown in Fig. 4). Consequently, the conformation e was selected as metastable conformation for the following experiments.

Close

Fig. 4

The RMSF profile of C-alpha in conformation b and conformation e. The RMSF profile of conformation b and conformation e was plotted in the upper and lower panel, respectively. The RMSF of C-alpha atom in A-loop was highlighted in red square.

Export to PPT

3.4 The identified metastable structural model showed enhanced efficacy targeting TKIs which was effective against D816V

A careful examination of the identified metastable conformation revealed it possesses the structural characteristics of both autoinhibitory c-Kit kinase and active c-Kit kinase. Although the A-loop adopted a partially extended conformation, the DFG motif preserved an in conformation, indicating it may represented a semi-active conformations, which may existed in the mutated c-Kit. Accordingly, a small dataset containing inhibitors targeting wild-type KIT (imatinib and masitinib) and KIT D186V (bezuclastinib, avapritinib, BLU-263, and elenestinib) was constructed, and their binding affinity of identified metastable conformation was predicted by molecular docking in standard precision mode. Interestingly, imatinib and masitinib, two inhibitors that were not effective on D816V mutations showed higher binding affinity towards conformation b than conformation e, while the selective KIT D816V inhibitors showed higher binding affinities towards conformation e(Fig. 5A). The detailed analysis of the binding pocket and binding poses revealed that the partial unfolding of the A-loop evoked the conformational alterations in the sidechain of V603, L595, F811, L799, and Y672, and led to the formation of a hydrophobic pocket adjacent to the DFG motif. At the same time, the rearrangement of the A-loop also evoked the conformational change of charged residues, such as D677, D816, K623, and R796. On the one hand, the newly formed hydrophobic pocket could stabilize the hydrophobic scaffold of inhibitors targeting KIT D816 mutations. On the other hand, the alteration of positively charged K263 and R796 could impair the repulse between KIT and inhibitors, and further stabilize the selective KIT D816V inhibitors within the binding pocket (Fig. 5 B and C).

Close

Fig. 5

The identified metastable structural model showed enhanced efficacy in targeting TKIs which was effective against D816V. The predicted binding free energy of identified metastable states and selected TKIs was obtained from molecular docking (Glide in standard precision mode). The docking score of each TKI was shown in A, and the three-dimensional structure of the complex was shown in B and C. The carbon atoms in the ligand were colored in white. The carbon atoms within conformation b and conformation e were colored cyan and pink, respectively. The nitrogen, oxygen, and chloride atom was colored blue, red, and green, respectively.

Export to PPT

3.5 Strucutre-based virtual screening identified vorolanib, semaxanib, henatinib, and pexmetinib may possess promising inhibitory effects on KIT D816V mutation

To identify the potential KIT D816V inhibitors from kinase inhibitors, a structure-based virtual screening process was performed, and for the compounds from PKIDB, a collection of kinase inhibitors was employed. The threshold of predicted binding free energy was set to −6.50 Kcal/mol, and we identified 114 inhibitors out of a total of 369 compounds. As shown in Supplementary Table 1, the majority of hit compounds belonged to the Tyr kinase family. Besides, we found our docking-based virtual screening protocol identified the majority of selective D816V inhibitors, such as toceranib, amuvatinib, avapritinib (predicted binding free energy: −7.537 Kcal/mol), and nilotinib (predicted binding free energy: 7.329 Kcal/mol) from available tyrosine kinase inhibitors. The ten top-scored compounds were listed in Table 1, and we found two out of ten selected compounds were selectively targeting D816V mutations. Consequently, we identified vorolanib, semaxanib, henatinib, and pexmetinib as potential inhibitors that target the D816V mutations.

3.6 In vitro validation of virtual screening

To evaluate the potential therapeutic effects of identified inhibitors, cell viability of the HMC1.1, and HMC1.2 cell lines, which possess D816V single mutation and D816V, V560G dual mutations was measured. As shown in Fig. 6, Vorolanib, the top-scored compound showed promising cytotoxicity against not only the HMC 1.1 cells, but also the HMC 1.2 cells. The results from the in vitro experiment further confirmed our prediction.

Close

Fig. 6

In vitro validation of virtual screening. The cytotoxic effects of effects of Vorolanib, the top-scored compound, on HMC1.1 and HMC1.2 cell line was measured by CCK-8 assay. The pooled data was visualized as mean ± S.E.M of four independent experiments. The one-way ANOVA and Tukey post-hoc test were employed for statistical analysis. *: <0.05, **: <0.01, ***: <0.001.

Export to PPT