Dasatinib-SIK2 Binding Elucidated by Homology Modeling, Molecular Docking, and Dynamics Simulations.

Salt-inducible kinases (SIKs) are calcium/calmodulin-dependent protein kinase (CAMK)-like (CAMKL) family members implicated in insulin signal transduction, metabolic regulation, inflammatory response, and other processes. Here, we focused on SIK2, which is a target of the Food and Drug Administration (FDA)-approved pan inhibitor N-(2-chloro-6-methylphenyl)-2-(6-(4-(2-hydroxyethyl)piperazin-1-yl)-2-methylpyrimidin-4-ylamino)thiazole-5-carboxamide (dasatinib), and constructed four representative SIK2 structures by homology modeling. We investigated the interactions between dasatinib and SIK2 via molecular docking, molecular dynamics simulation, and binding free energy calculation and found that dasatinib showed strong binding affinity for SIK2. Binding free energy calculations suggested that the modification of various dasatinib regions may provide useful information for drug design and to guide the discovery of novel dasatinib-based SIK2 inhibitors.

Introduction

Salt-inducible kinases (SIKs) control cyclic adenosine monophosphate (cAMP)-dependent production of the anti-inflammatory cytokine interleukin-10 (IL-10) in macrophages.[1−6] The three SIK isoforms reported are SIK1, SIK2 (QIK), and SIK3 (QSK).[7] SIK2 modulates metabolic pathways,[2] steroidogenesis,[8] adipogenesis,[7] adipocyte energy metabolism,[9] fatty acid oxidation,[10] and centrosome splitting.[11] SIK2 overexpression has been implicated in the development of chronic inflammatory diseases,[12] gastric cancer,[13,14] and acute kidney injury.[15] Hence, it is an appealing pharmacological target for the treatment of macrophage-driven diseases.[10,11,13,16−19]

Recently, numerous efforts have focused on the development of SIK modulators, especially those that target SIK2.[13,20−27] SIK inhibitors under development include HG-9-91-01, YKL-06-062, MRT67307, MRT199665, and ARN-3236 (Figure ). HG-9-91-01 occupies adenosine triphosphate (ATP)-binding sites and small hydrophobic vesicles near them.[28] YKL-06-062 is a selective SIK inhibitor.[29] ARN-3236 was the first orally administered SIK2 inhibitor.[30] High-throughput screening assays and probes have been developed to identify SIK inhibitors.[24,31] Thus far, however, no clinically approved inhibitors specifically targeting SIKs have been reported. Certain Food and Drug Administration (FDA)-approved kinase inhibitors such as crizotinib, N-(2-chloro-6-methylphenyl)-2-(6-(4-(2-hydroxyethyl)piperazin-1-yl)-2-methylpyrimidin-4-ylamino)thiazole-5-carboxamide (dasatinib), erlotinib, gefitinib, and pazopanib inhibit SIK functions.[13,32,33] Nevertheless, they target several other protein kinases as well. Therefore, the development of novel SIK2 inhibitors is necessary.

Figure 1

Structures of SIK2 inhibitors. (A) Structures of dasatinib, ARN3236, YKL-06-062, HG-9-91-01, MRT67307, and MRT199665. (B) Labeling of the SIK2 structure.

The lack of small molecule inhibitors targeting SIK2 emphasizes the necessity of exploration and identification of novel drug candidates. Herein, we describe the mechanism by which the FDA-approved kinase inhibitor dasatinib binds SIK2. N-(2-Chloro-6-methylphenyl)-2-(6-(4-(2-hydroxyethyl)piperazin-1-yl)-2-methylpyrimidin-4-ylamino)thiazole-5-carboxamide (dasatinib; BMS-354825; Sprycel; Figure ) is a first-line drug used for the treatment of BCR-ABL-positive leukemias.[34] It can also be used to treat chronic myelogenous leukemia (CML), lymphomas, and advanced prostate and breast cancers.[14,35−40] Adverse events, such as gastrointestinal disorders, hemorrhage, and endothelial permeabilization, leading to the development of peripheral edema and pleural effusion were reported in certain patients who were administered with dasatinib.[37,41−43] Dasatinib is a broad-spectrum kinase inhibitor targeting certain protein kinases such as YES proto-oncogene 1 (YES1), proto-oncogene tyrosine-protein kinase SRC (SRC), tyrosine-protein kinase LYN (LYN), FYN related Src family tyrosine kinase (FRK), FYN, EPH receptor B2 (EPHB2), EPH receptor A2 (EPHA2), discoidin domain receptor tyrosine kinase 1 (DDR1), Abelson leukemia virus tyrosine kinase (ABL)2, receptor interacting protein kinase (RIPK)2, LIM kinase 1 (LIMK1), and SIK2.[14,44,45] Dasatinib noncovalently binds and inhibits SIK2 functions.[14] In vitro, it inhibits functions of SIK isoforms with half-maximal inhibitory concentration (IC50) in a nanomolar range (<3 nM for SIK1, <3 nM for SIK2, and 18 nM for SIK3).[27,46] As it exhibits remarkable activity against SIKs and is also a pan inhibitor against kinases, dasatinib may be considered a lead compound in the design and/or development of new, highly selective SIK2 inhibitors.

A reliable receptor-ligand structure is vital for knowledge- or structure-based drug development. Elucidation of drug–target binding requires an understanding of ligand-receptor incorporation. Homology modeling (HM) is a widely accepted computational method used for the prediction of protein structures.[47,48] It can clarify protein three-dimensional (3D) structures based on amino acid sequences.[49,50] Thus, it helps identifying novel drug candidates and may enable drug discovery in a faster, easier, cheaper, and more practical manner.[47,48,51−53] The new HM technique is widely used in drug development.[54,55] Molecular docking creates a static image of the drug–target complex and is also used in drug design.[56−64] Molecular dynamics (MD) helps elucidate the interactions in drug binding.[65−73]

Presently, no crystalline structural data have been reported for SIK2. This information gap hinders the development and improvement of SIK2 modulators. In the present study, SIK2 was constructed by performing HM, a dasatinib–SIK2 binding model was constructed by molecular docking, and extensive MD simulations were performed. Binding free energies were calculated using the molecular mechanics generalized Born (GB) surface area (MM/GBSA) method, which has proved to be a powerful and valuable tool.[74−76] Certain residues involved in dasatinib–SIK2 binding were determined by analysis of per-residue energy decomposition. We believe that these simulations may provide useful information to facilitate the design of innovative SIK2 inhibitors.

Results and Discussion

Three-Dimensional SIK2 Structure

We used HM to plot a 3D structure for the SIK2 kinase domain (KD). We used the SIK2 sequence (residues 20–271) as the query option and found that data on >100 crystal structures extracted from the PDB showed high homology with SIK2. The BLASTP criteria comprised E value <10–50 and query cover >99%. SIK2 and microtubule affinity-regulating kinase (MARK) are members of the calcium/calmodulin-dependent protein kinase (CAMK)-like (CAMKL) family. Hence, high kinase domain (KD) homology is expected between SIK2 and MARK. MARK1–MARK4 were selected to construct the 3D SIK2 structure. We selected only chain A of the four crystalline structures with PDB IDs 2HAK(77) (MARK1), 5EAK(78) (MARK2), 2QNJ(79) (MARK3), and 5ES1(80) (MARK4) and used them as templates (Figure S1). The similarity between SIK2 and MARK1–4 was >62%. Four high-homology templates were used to construct the SIK2 structures.

Four 3D SIK2 protein structures were plotted using SWISS-MODEL according to the four MARK templates. The structures are depicted in the superimposition in Figure and labeled as SIK2-I (MARK1), SIK2-II (MARK2), SIK2-III (MARK3), and SIK2-IV (MARK4). Figure shows that most of the regions in the four structures presented with good overlap. A ProCheck Ramachandran plot was used to assess the structure reliability (Figure S2). Fewer than two residues were located in the disallowed regions. The proportions of residues in most of the favored regions were 89.7, 90.6, 91.5, and 90.1% for SIK2-I, SIK2-II, SIK2-III, and SIK2-IV, respectively. The global model quality estimation (GMQE) scores were 0.79, 0.78, 0.78, and 0.77 for SIK2-I, SIK2-II, SIK2-III, and SIK2-IV, respectively. However, QMEAN results revealed that certain low-score regions (<0.6) were positioned for these models (Figure S3). In the MARK templates, the low-score regions represented various conformations and were relatively more flexible. The foregoing scores were generally better than the commonly accepted standard for HM.[47,48,51,81−83] Therefore, the constructs served as templates in subsequent simulations.

Figure 2

Diagram of the four SIK2 models. (A) Overall structure. (B) T-loop highlight.

Molecular Docking

Before performing molecular docking to construct the initial dasatinib/SIK2 binding model, we evaluated the docking power via a “re-docking” strategy. We selected a previously reported crystal structure for the ABL kinase domain (PDB entry code 2GQG(84)), which formed complexes with dasatinib. The root-mean-square deviation (RMSD) obtained between the crystal structure and the conformer with the lowest energy (−11.92 kcal/mol) based on docking results was 1.0 Å (Figure S4). Hence, we performed molecular docking to construct the dasatinib/SIK2 complex.

A PDB survey revealed that ≥20 kinase crystal structures formed complexes with dasatinib (Table S1).[84−97] All structures demonstrated nearly the same binding pattern with dasatinib, as highly conserved binding pockets were present for all of the structures. Figure S5 shows similar geometric dasatinib conformers binding with various protein kinases. Dasatinib could form a close binding interaction with SIK2. We plotted four molecular docking models with the lowest binding energies (Figure ). The ABL crystal structure was included for comparison. A relatively good overlap was observed. Therefore, the binding pockets observed in the SIK2 homology model may be suitable for consideration and the status of dasatinib in the binding sites should be reliable. Figure S6 shows the presence of putative interactions between dasatinib and the SIK2-I model. A hydrogen bond network that potentially facilitated the binding affinity was observed between dasatinib and SIK2. Certain hydrogen bonds are plotted in Figure S5A, including those between dasatinib N1 and T96 and between A99, N2, and N3. The atom labels for dasatinib and the T96 and A99 residues of SIK2 are illustrated in Scheme S1. The foregoing analyses are exclusively based on the docking structures obtained herein. Long-term molecular dynamics simulations should be performed to establish system stability and ligand recognition.

Figure 3

Selected results of dasatinib docking with SIK2. (A) SIK2-I, (B) SIK2-II, (C) SIK2-III, and (D) SIK2-IV. The ABL crystal structure (PDB ID: 2GQG and color with yellow) was included for comparison.

Molecular Dynamics Simulation

The RMSD values for the SIK2-I backbone atoms were calculated according to the structures obtained by docking and have been illustrated in Figure A. The RMSD was 1.91 ± 0.17 Å throughout the entire 500-ns MD simulation. Small final fluctuations indicated that the entire system was stabilized. The ligand remained securely in the pocket. Similar results were obtained for the other three models (Figure S7). All four HM structures were sufficiently stable for subsequent analyses. The root-mean-square fluctuation (RMSF) values were analyzed based on the 500-ns MD trajectory information. RMSF is based on the fluctuation of all amino acids and quantifies the stability of specific residues during MD simulation. The amino acid residues at positions 39–43, 54–59, 80–85, 116–120, 188–194, and 218–222 and the C-terminal 254–258 of SIK2-I presented with large fluctuations (RMSF > 1.00) (Figure B). The same phenomena were observed for the other three models (Figure S8). The relatively flexible T-loop conformations were in good agreement with HM and other kinases such as MARK2[98] and MARK4.[99] The T-loops of RIPK family members may exhibit open or closed conformations at the active site and can recognize various inhibitors.[100] Here, we identified open and closed conformations in our HM-based constructs. They were indicated by the snapshots obtained at 100, 200, 300, 400, and 500 ns from the MD trajectory and have been aligned and depicted for dasatinib/SIK2-I in Figure C. The other models are illustrated in Figures S9 and S10. Only the piperazine group position on the ligand showed fluctuations (Figures C and S9). This generic fragment conformation in the solvent environment can improve the bioavailability.[101−103] Therefore, the SIK2 conformation should be stable when dasatinib is at the ATP-binding site.

Figure 4

(A) Plot of RMSD of Cα atom vs time for the 500-ns MD simulation of SIK2-I. (B) RMSF variations of the Cα atom of SIK2 from the 500-ns MD simulation of SIK2-I. (C) SIK2-I snapshots along the dynamic simulation timelines at 100, 200, 300, 400, and 500 ns. For clarity, the water molecules have been removed. Dasatinib is represented by a stick diagram, while SIK2 is depicted by a cartoon diagram.

The initial docking models identified a substantial hydrogen bond network present between the receptor and the ligand. In contrast, dasatinib possesses active hydrogen bond acceptors and donors. The hydrogen bonds formed between dasatinib and SIK2 were then examined. Hydrogen bond occupancy analyses over the MD trajectories for all four binding complexes are summarized in Figure . Only a few hydrogen bonds formed between the ligand and the protein. However, certain hydrogen bond occupancies >86% and the hydrogen bond number were also detected (Figures and S11). Dasatinib N1 formed hydrogen bonds with T96 OG1. The latter is a gatekeeper residue for protein kinases. The reduction in the inhibition of the HG-9-91-01 strain against T96Q was >500 times greater than that for the wild-type SIK2.[28,104] Two other hydrogen bonds were observed in the dasatinib hinge loop and were also found for other protein kinases.[85,86,105]

Figure 5

(A) Hydrogen bond network analysis of interactions between dasatinib and SIK2 for SIK2-I, SIK2-II, SIK2-III, and SIK2-IV. Occupancy is expressed as percent of the period (500 ns) during which specific hydrogen bonds are formed. (B) Hydrogen bond network analysis of the interactions between dasatinib and SIK2-I. (C) Statistical hydrogen bond number profile along the 500-ns MD simulation for SIK2-I. Hydrogen bond is defined as the distance between the acceptor and donor atoms <3.5 Å, with an internal angle between the H-acceptor and H-donor >120°.

Binding Free Energies

Binding models generated using the docking algorithm and MD simulations depict the interactions between dasatinib and the SIK2 protein. However, they cannot be used to identify the best model. Here, we used MM/GBSA to calculate the absolute dasatinib–SIK2 binding free energies for all four systems (Table ). Uncertainties shown in parentheses were calculated as the root-mean-square error for all frames extracted from the trajectories. All TΔS and enthalpies were negative (<−24.00 and <−46.00 kcal/mol, respectively). Thus, binding complex formation is an enthalpy-driven process. For SIK2-I, the calculated binding free energy was ∼−29.37 kcal/mol. For SIK2-II, SIK2-III, and SIK2-IV, the ΔGbindcal values were −20.09, −21.97, and −22.70 kcal/mol, respectively. Therefore, SIK2-I was probably the preferred model as its T-loop pointed toward the ATP-binding site, which represents a closed conformation for the T-loop. In other words, the binding free energies of dasatinib to the SIK2 largely depend on the conformation of the T-loop. Specifically, in the model of SIK2-I, dasatinib formed a hydrogen bond between the O1 atom of dasatinib and the side chain of S169 in the T-loop for occupancy with 61.60%. This hydrogen bond was not found in the other three models. However, we could not exclude the existence of the other three models as the kinase conformation could be changed to fit and bind the inhibitor. For example, dasatinib can bind active and inactive ABL1 conformations and form a strong binding with ABL1.[106−108]

Table 1

Binding Free Energies (ΔGbindcal) for Dasatinib/SIK2 Complexes and Decomposition of Electrostatic Interaction (Eele), van der Waals (vdW) Interaction (EvdW), Solvation Free Energies (EGB), and Entropy (TStotal)a

	SIK2-I	SIK2-II	SIK2-III	SIK2-IV
EvdW	–60.14 (2.87)	–54.21(2.87)	–53.56(2.91)	–54.63(3.15)
Eele	–25.21(6.14)	–26.08(5.60)	–24.31(5.35)	–24.02(5.98)
EGB	36.99(4.58)	38.81(5.21)	37.03(4.97)	35.87(5.43)
Esurf	–6.63(0.23)	–6.00(0.28)	–5.93(0.29)	–5.97(0.31)
Ggas	–85.35(6.32)	–80.29(6.26)	–77.88(6.15)	–78.65(6.67)
Gsolv	30.36(4.56)	32.81(5.14)	31.09(4.89)	29.91(5.34)
Egas + Gsol	–54.99(3.39)	–47.48(2.98)	–46.78(3.15)	–48.74(3.12)
TΔStotal	–25.62(5.22)	–27.39(5.29)	–24.81(4.76)	–26.04(5.50)
ΔGbindcal	–29.37(6.22)	–20.09(2.98)	–21.97(5.71)	–22.70(6.32)

Energy values are presented in kcal/mol. Uncertainties shown in parentheses were calculated as the root-mean-square error for all frames extracted from the trajectories.

The polar (Eele + EGB) and nonpolar (EvdW + Esurf) terms usually affect inhibitor binding. The polar terms were ∼11.78, ∼12.73, ∼12.72, and ∼11.85 kcal/mol for SIK2-I, SIK2-II, SIK2-III, and SIK2-IV, respectively. Electrostatic interactions are nullified by the solvation effect. Positive values for the electrostatic contribution indicate that polar interactions between the ligand and the receptor are antagonistic to this binding. In contrast, van der Waals interactions are conducive to binding as they constitute the main nonpolar contribution. The total nonpolar values were −66.77, −60.21, −59.49, and −60.60 kcal/mol for SIK2-I, SIK2-II, SIK2-III, and SIK2-IV, respectively. In simulations conducted in the present study, the nonpolar values were markedly greater for SIK2-I than the others. This finding was consistent with our conclusion that SIK2-I exhibited a stronger binding affinity than SIK2-II, SIK2-III, or SIK2-IV. Hence, hydrophobic interaction is the predominant factor responsible for dasatinib–SIK2 binding.

Free Energy Decomposition

The calculated binding free energy showed that the nonpolar term played the most important role in complex formation. The per-residue free energy decomposition strategy enables the analysis of the inhibitor–protein interaction. The interaction energy between each residue and the inhibitor was computed using the MM/GBSA decomposition protocol.

Certain hydrophobic residues possessed substantial subtotal binding free energies. For SIK2-I, L26 in the P-loop demonstrated the most favorable contribution (Table ). P-loop contribution was also detected for TYK2 with its inhibitors[109] and imatinib with ABL1.[110,111] T96 and A99 each contributed <−2.00 kcal/mol because of the hydrogen bonds formed with dasatinib. Similar findings were obtained for HG-9-91-01 and KIN112 binding with SIK2.[24,104,112−114] Previous empirical and theoretical studies on protein kinases reported observations in similarity to those for V981 in the TYK2 hinge loop.[109] L26, T96, and A99 are the most favorable residues for SIK2-II, SIK2-III, and SIK2-IV binding and are, therefore, vital to dasatinib binding. V34, A47, K49, Y98, G102, and L149 each contributed <−1 kcal/mol. Moreover, the nonpolar terms and, especially, the van der Waals interactions were major components of the subtotal binding free energies for most residues (Tables S2–S5). The ΔEele term of K49 was ∼−5.23 kcal/mol. However, the electrostatic term (ΔEele + ΔGsol,GB) was <0.50 kcal/mol. Thus, the solvent effect counteracted the amino contribution for SIK2-I (Table S2). The electrostatic effect of S169 (−0.56 kcal/mol) nearly equaled the nonpolar contribution of −0.47 kcal/mol (ΔEvdW + ΔGsol,np). D160 with a negative charge contributed 0.79 kcal/mol to SIK2-I and 0.60, −0.25, and −0.05 kcal/mol to SIK2-II, SIK2-III, and SIK2-IV, respectively. Hence, polar interactions markedly contributed to the binding of certain residues even though the overall electrostatic interactions were unfavorable.

Table 2

Free Energy Decomposition of the Dasatinib/SIK2 Complexes at the Level of the Contributions of Individual Residuesa

residues	SIK2-I	SIK2-II	SIK2-III	SIK2-IV
L26	–2.75	–2.38	–2.37	–2.37
V34	–1.74	–1.49	–1.48	–1.57
A47	–1.06	–1.08	–1.13	–1.07
I48	–0.69	–0.66	–0.65	–0.67
K49	–1.84	–0.63	–0.53	–0.63
M71	–0.49	–0.43	–0.59	–0.64
I80	–1.33	–1.18	–1.18	–0.80
L94	–0.65	–0.67	–0.71	–0.60
V95	–0.74	–0.73	–0.58	–0.86
T96	–2.20	–2.44	–2.10	–2.78
E97	–0.73	–0.86	–0.84	–0.77
Y98	–1.57	–1.68	–1.78	–1.64
A99	–2.14	–2.37	–2.01	–2.43
N101	–0.66	–0.84	–0.75	–0.83
G102	–1.44	–1.45	–1.52	–1.39
E103	–0.29	0.35	0.28	0.13
D106	0.19	0.27	0.25	0.24
L149	–1.86	–1.71	–1.84	–1.79
D160	0.79	0.60	–0.25	–0.05
S169	–1.02	0.02	0.02	0.00

Energy values are presented in kcal/mol.

The main contributions of T96 and A99 to SIK2-I, SIK2-II, SIK2-III, and SIK2-IV are attributed to the hydrogen bond network formed between the inhibitor and the protein. Residues responsible for significant contributions to overall binding were situated mainly at the inhibitor-binding groove. Residues forming hydrogen bonds with the dasatinib core were conserved in the protein kinase.[84−97] Therefore, the hydrogen bonding network could be vital for ligand binding but might not be responsible for selectivity. L26, V34, A47, K49, Y98, G102, and L149 contributed mainly to dasatinib–SIK2 binding. However, these residues could fit the inhibitor conformations to increase selectivity. The influence of van der Waals interactions on hydrophobic amino acid binding may be essential for dasatinib–SIK2 binding.

Simulated Annealing Molecular Dynamics (SAMD)

Kinase activation loops have two conformers that differ in terms of their ability to recognize various inhibitors. Imatinib could bind only the inactive ABL1 conformation, whereas dasatinib could also bind the active ABL1 conformation.[84,106−108,115,116] Therefore, dasatinib targets imatinib-resistant ABL and increases the binding affinity. HM and MD revealed the closed (SIK2-I) and open (SIK2-II, SIK2-III, and SIK2-IV) conformations of T-loop. Dasatinib bound SIK2 in all four simulation models. However, classical MD simulations and binding free energy calculations exhibited different binding affinities. Therefore, additional calculations should be performed to clarify the manner in which the T-loop influences the inhibitor–SIK2 interactions. As large-scale conformational changes occur, classical MD with limited simulation time cannot be used to map the entire process. For this reason, a simulated annealing MD approach was applied. SAMD was used to successfully elucidate the protein–ligand complex binding model and conformational changes in proteins.[117,118]

We used the SAMD protocol to generate various binding conformations for the dasatinib/SIK2 complex. RMSD, temperature, and energy variation for the first SIK2-III simulation cycle are illustrated in Figure S12. During annealing, the system overcomes energy barriers and reaches different local minima on the potential energy surface. After SAMD was run for all four models, 400 conformations were clustered via hierarchical agglomeration and structural similarity patterns were identified (Figure A). Nine clusters represented the 400 conformations. The top four clusters accounted for >89% of all conformations. There were 100 (25%), 100 (25%), 100 (25%), and 59 (14.7%) conformations for SAMD-I, SAMD-II, SAMD-III, and SAMD-IV, respectively. The representative frames for the top four clusters are shown in Figure C–F. The T-loop was closed for SAMD-I, SAMD-II, and SAMD-IV. The distances between the T175 Cα atom and the dasatinib pyrimidine were 15.8, 13.6, and 11.5 Å for SAMD-I, SAMD-II, and SAMD-IV, respectively. However, the distance was 24.9 Å for SAMD-III. Hence, it should be labeled as an open conformer. Our SAMD results clearly show that when dasatinib is bound within the SIK2 active site, both close and open conformations can be found for the T-loop. In other words, the inhibitory activity of dasatinib against SIK2 might be maintained at both closed and open conformations of the T-loop. Of course, the closed T-loop can block the ATP access channel because the ATP-binding cleft is occupied by dasatinib and the T-loop.[112,119]

Figure 6

Cluster analysis of the 400 conformations for SAMD (A). T-loop conformation aligned for the representative structures of the top four clusters (B). Conformations for SAMD-I (C), II (D), III (E), and IV (F). Dasatinib is depicted using a stick diagram, while SIK2 is represented by a cartoon diagram.

Atomistic details of the contacts and noncovalent dasatinib–SIK2 interactions from SAMD are illustrated in Figure S13. Dasatinib forms several hydrogen bonds with SIK2 in four representative conformers extracted from SAMD (Figure S14). The energies of these hydrogen bonds were <−2.00 kcal/mol. A99 N formed hydrogen bonds with dasatinib N2 and contributed only −2.1 kcal/mol for SAMD-III but −3.9, −4.7, and −4.9 kcal/mol for SAMD-I, SAMD-II, and SAMD-IV, respectively. This hydrogen bond formation weakens during the transformation of closed to open conformers. The hydrogen bonds formed with A99 were generally stronger than those formed with T96 (Figure S14). However, hydrophobic interactions with L26, V34, A47, M71, I80, L94, A99, L149, and A159 were observed in all four cluster models. I48, L82, V95, and A174 established interactions with dasatinib in only one SAMD model. The free energy decomposition calculations suggested that these were key residues.

The foregoing results are consistent with the MD simulations and binding free energy calculations. The overall dasatinib/SIK2 interaction patterns were similar for all models. Dasatinib/SIK2 recognition should be dominated by hydrophobic interactions and hydrogen bonds. T-loop fluctuations cannot markedly alter the binding model.

Design Strategies

We can divide the dasatinib molecule into four regions based on the results obtained from the MD simulations and binding free energy calculations (Figure ). The 2-chloro-6-methylphenyl ring probes the hydrophobic pocket and is labeled as the R1 region. The dasatinib pyrimidine ring is near the P-loop and is designated as the R2 region. The piperazine group is identified as the R3 region pointing toward the solution. The dasatinib aminothiazole ring forms two strong hydrogen bonds with the SIK2 hinge loop and is defined as the R4 region. Seven new molecules were designed by modifying these four regions (Figure ). In particular, H1 was constructed without the methyl group in the pyrimidine. The H2 molecule was generated by removing the aminothiazole ring. Piperazine was deleted to generate the S1 compound. G1–G4 were obtained using different substitution groups to replace the 2-chloro-6-methylphenyl group. Docking, MD simulations, and binding free energy calculations were performed to establish the binding affinities for the novel molecules. This approach may elucidate the contributions of the various regions of dasatinib to its overall inhibition capacity.

Figure 7

Structures of designed analogues and binding free energy of SIK2 based on MM/GBSA calculations with final 50-ns MD simulations.

We adopted the SAMD-I protein structure as the receptor for the newly designed compounds to construct complex structures by docking. We used previously described computational protocols. The complex models are illustrated in Figure S15. First, because the H2 molecule lacks an aminothiazole ring as shown in Figure , it would further block interactions with the hinge loop and thus impair the ligand binding. Earlier experimental studies corroborated this finding for ABL kinases.[34,84−86] Therefore, we did not run the simulation to evaluate the H2 molecule. G2, G3, and G4 were similar in terms of ligand binding. Since G3 can form the most stable (−10.48 kcal/mol) complex with SIK2 than the other two systems (−9.80 and −9.60 kcal/mol for G2 and G4, respectively), we then simply carried out simulations only for G3 rather than G2 and G4. Finally, 100-ns MD simulations were run for the SIK2 system complexed with S1, H1, G1, and G3, respectively. RMSD (Figure S16) and binding free energies (Figure ) were also calculated. The systems were generally stable after 50 ns.

Close inspection of the MD trajectories for all four complexes revealed that S1, H1, and G3 were tightly bound (Figures S17–S19). Nevertheless, the entire molecule moved toward the solution phase for G1/SIK2 binding (Figure S20). In G1, the 2-chloro-6-methylphenyl ring in R1 was replaced by a cyclopropane. The simulation performed in the present study lowered the binding affinity. It demonstrated a binding free energy of −9.22 kcal/mol as opposed to a range of −29.37 to −20.09 kcal/mol for dasatinib/SIK2. The dasatinib aromatic ring may not be replaced by cyclopropane. The H1/SIK2 binding free energy was −26.08 kcal/mol, which was ∼3.3 kcal/mol higher than that of the dasatinib/SIK2-I complex. Nevertheless, the former was superior to that of the other dasatinib/SIK2 models. The entire H1 molecule bound well in the ATP-binding site of SIK2. ΔG(S1/SIK2) was −20.29 kcal/mol, which was ∼9.08 kcal/mol higher than that of the SIK2-I model. However, ΔG(S1/SIK2) was probably comparable to those of SIK2-II, SIK2-III, and SIK2-IV. Previous studies showed that the dasatinib piperazine fragment exerted negligible influence on inhibitor binding.[101,103] Tiwari et al. replaced the terminal hydroxyethyl group with amino and fatty acids to enhance bioavailability.[101] Our simulations provided theoretical insights into dasatinib modification. For the G3 molecule, the binding free energy was −22.58 kcal/mol. Binding affinity was not markedly affected when dasatinib established interactions with the SIK2 hydrophobic pocket modified by a nonaromatic ring.

These four small molecules, S1, H1, G1, and G3, were designed based on the molecular framework of dasatinib, which could shed light on the binding affinity of different parts of dasatinib with SIK2. In this way, such kind of model study could provide some useful information for future inhibitor designs based on dasatinib. Dasatinib molecules with substitutions at S1, H1, and G3 exhibit a remarkable binding affinity for SIK2. The modeling and structuring strategies used herein helped identifying potentially optimal dasatinib modifications. Replacement of R1 with a six-membered ring, R2 with an aromatic ring, and R3 with a polar group can enhance the bioavailability and binding affinity between the ligand and SIK2.

Computational Details

Homology Modeling

Human SIK2 protein sequence (UniProt ID: Q9H0K1) was extracted from the UniProt database.[120] The full SIK2 protein consisted of 926 amino acids and was composed of a kinase domain near the N-terminal (KD), a central sucrose nonfermenting (SNF) protein kinase homology domain (SNH), and a phosphorylation domain near the C-terminal (SIKC).[1,2,15] KD was the substrate-binding site during SIK2 phosphorylation.[6] Most inhibitors interacted with the protein KD.[23,24,29,30,121] Here, only KD (residues 20–271) was selected for model construction. A BLASTP[122,123] search against the protein data bank (PDB; http://www.rcsb.org/pdb/)[124] was conducted to determine homologous proteins. Four protein crystal structures for microtubule affinity-regulating kinase 1 (MARK1) (PDB ID: 2HAK_A[77]), MARK2 (PDB ID: 5EAK_A[78]), MARK3 (PDB ID: 2QNJ_A[79]), and MARK4 (PDB ID: 5ES1_A[80]) were selected for SIK2 homology modeling. These proteins are members of the CAMKL family and share high sequence similarity.[125,126] However, the T-loop region in these four crystal structures has different conformations. At the same time, there are some missing residues among the T-loops of these four crystal structures. Therefore, the inclusion of these four proteins could provide more structural information from our homology modeling. The online SWISS-MODEL service tool[127,128] was used to construct the 3D SIK2 structures. The SIK2 homology modeling structures were evaluated by GMQE[127] and a Ramachandran plot via ProCheck.[129,130]

Molecular Docking

The receptor–ligand complex structure was constructed by performing molecular docking according to the data on the SIK2 receptor protein obtained by homology modeling. The protein structure was pretreated using AutoDockTools v.1.5.6[131] and included hydrogen addition, Gasteiger charging,[132] unreasonable atomic overlap adjustment, and so on. The dasatinib structure was optimized via the PM3 force field[133,134] by MOPAC2016 (James J. P. Stewart, Stewart Computational Chemistry, Colorado Springs, CO; http://OpenMOPAC.net (2016)). A 30 × 60 × 30 grid map with 0.375 Å grid spacing was generated using AutoGrid v.4.2[135] and was based on the center of the active SIK2 groove. One hundred conformations per system were generated using the Lamarckian genetic algorithm in AutoDock.[135] Optimal conformations were acquired for all docking models.

Molecular Dynamics Simulation

The dasatinib/SIK2 complexes were obtained by performing molecular docking. General Amber force field (GAFF)[136] generation was considered to describe the dasatinib force field. The geometric structure was optimized at the B3LYP/6-31G theory level using Gaussian 09.[137] To obtain data on the partial atomic charges, the restrained electrostatic potential (RESP) protocol[138] was implemented at the B3LYP/6-31G* theory level. The force field parameters were generated using the Antechamber module. Further, the AMBER ff14SB force field[139] was used to create the protein topology parameters. The complexes were dissolved in TIP3P water[140] in silico box with ∼93 Å × ∼85 Å × ∼87 Å. The systems were neutralized using chloride (Cl–) ions. The final system included 4195 solute and ∼17 000 solvent (water) molecules. Periodic boundary conditions were set to avoid edge effects, and a 12-Å radius cutoff was established to enable van der Waals interactions. The particle mesh Ewald (PME) algorithm[141] was used to calculate long-range electrostatic interactions. The SHAKE algorithm[142] was used to constrain the covalent bonds associated with the hydrogen atoms. To reduce unfavorable interactions caused by solvents and ions, the system was subjected to 9000 steepest descent method steps followed by a 1000-step conjugate gradient. All solute molecules were fixed at the initial position. Thereafter, a 10 000-step conjugate gradient was used to optimize the entire system including the solute and solvent molecules. After the performance of the initial two-step system minimization, the overall system temperature was increased from 0 to 300 K in 200 ps via Langevin dynamics and collision frequency g = 2.0 ps–1. Pressure was maintained at 1 bar for 200 ps using isotropic position scaling[143,144] and temperature = 300 K. The system was then equilibrated at 300 K and 1 bar over 20 ps within the NPT ensemble. The entire system was subjected to 500-ns MD simulations for data collection and analysis. The integration step size was set to 2 fs throughout the MD process. All dynamics were performed using the CUDA version of PMEMD in AMBER 12.[145] The CPPTRAJ module[146−148] was used to analyze the MD trajectory data.

Binding Free Energy

Molecular mechanics generalized Born surface area (MM/GBSA)[149−151] was used to calculate binding free energies to quantitate the binding affinity of dasatinib for SIK2. MM/GBSA is represented by the following equationswhere ΔGcomplex, ΔGprotein, and ΔGligand indicate the free energies of dasatinib/SIK2 as a complex, SIK2 as a protein, and dasatinib as a ligand. MM/GBSA was used to calculate these values based on the statistical average obtained for one MD trajectory. The solvation free energy (Gsol) comprises electrostatic (Gel) and nonelectrostatic (Gnonel) terms. Gnonel represents the combined effects of the unfavorable cost of surface formation and favorable van der Waals (vdW) interactions between the solute and solvent. This factor was evaluated using the equation γ = SA + b, where γ = 0.0072 kcal/Å2 and b = 0.0 kcal/mol. The solvent-accessible surface area (SA) was estimated by linear combinations of pairwise overlap (LCPO).[152]Gel was calculated using the generalized Born (GB) equation.[153,154] The interior of the system was assumed to be filled with dielectric constant (ε) of 1. The solvent water considered had a high dielectric constant (ε = 80). The gas-phase energy contributions (Egas) were calculated using mmpbsa_py_energy in the AmberTools package according to the force field with mbondi2 for PBRadii.[155] For each complex system, the binding energies were averaged over 1000 frames of the last 200-ns MD trajectory. Entropy contributions to the binding free energy were added as a further refinement and averaged for a 2-ns interval of the last 200 ns of MD. The entropy was estimated based on the changes in the degrees of freedom including translation, rotation, and vibration.[156] The free energy was calculated using MMPBA.py.[157]

The contribution of each residue to the binding free energy was evaluated at the atomic level by considering free energy decomposition.[158] The per-atom contribution was summed over the residue atoms, such that the contributions of the individual residues could be “non-mutably” estimated at the atomistic level. The same generalized Born model was used to determine the electrostatic contribution to the solvation energy. Internal energy calculation was excluded under the assumption of a single trajectory simulation. The corresponding nonpolar solvation energy per atom was obtained based on the corresponding SA. No entropy contribution was included in this approach. The total relative binding free energy of a given residue was obtained by summing the contribution of each atom present in it. The separate backbone and side-chain contributions were determined based on the analysis of the relevant atoms. MM/GBSA free energy decomposition was used to consider the same snapshots extracted for the binding free energy calculation.

Simulated Annealing Molecular Dynamics

The activation (T-) loop plays an important role in kinase activity.[159−163] Simulated annealing molecular dynamics (SAMD) was conducted to determine the most probable T-loop conformation when dasatinib bound SIK2. The final MD simulation snapshots for the complex systems were retained as initial structures. The dasatinib/SIK2 SAMD simulations were only performed in vacuo when all water molecules and chloride ions were removed. The system was minimized and equilibrated using the aforementioned computational protocol. After equilibration, 100 SAMD cycles were run per system. The heavy backbone atoms were restrained with 2.0 kcal/(mol·Å2) to maintain the equilibrated conformation. Meanwhile, the T-loop atoms freely moved during the simulation. Each cycle was commenced from the same conformation at randomized velocities to enhance the conformational space diversity. The system temperature was linearly increased from 300 to 1200 K over 10 ps during the heating phase. The system was then relaxed for ∼15 ps at 1200 K and cooled to 300 K within 25 ps. Another 25-ps relaxation dynamic simulation was run at 300 K. The final conformations of each cycle were considered for data analysis. Four hundred frames were obtained for the four systems; they were analyzed by hierarchical agglomeration and visualized using PyMOL 2.1[164] and VMD 1.8.[165]

Conclusions

Therapeutic agents with target protein specificity are highly attractive in drug design. However, the lack of detailed information on receptor–ligand interactions restricts development in this field. Crystallographic structural data are invaluable in pharmaceutical research. For certain molecules such as SIK2 investigated here, data on only amino acid sequences are available in the literature. Systematic calculations based on molecular docking, molecular dynamics, and binding free energy calculations should be conducted to develop specific inhibitors for pharmaceutical targets such as SIK2.

In the present study, the SIK2 3D structure was constructed using high-resolution MARK family templates and the FDA-approved kinase inhibitor dasatinib was used as the template ligand. Dasatinib–SIK2 interactions were investigated by performing molecular docking, MD simulation, and binding free energy calculations. Certain compounds were designed to quantify the contributions of various regions of dasatinib to its SIK2 binding affinity. Our simulation showed that dasatinib exhibited a linear structure in the SIK2 binding pocket, as observed in studies on other protein kinases. Hence, dasatinib nonspecifically and directly inactivates the protein kinase. In contrast, the ligand occupies the purine binding site with a hinge loop and extends to the bottom of the pocket by interacting with the hydrophobic wall formed by I80, L94, A159, and M71. The T-loop alters the conformation to fit the inhibitor. Thus, inhibitors may indirectly influence protein–protein interactions between SIK2 and its substrates. The simulations in the present study demonstrate that A99 and T96 form hydrogen bonds with the aminothiazole and amide fragments in the ligand, which, in turn, fix the conformation in the ATP-binding site. Free energy decomposition and simulations based on annealing MD indicated that L26, V34, A47, M71, I80, L94, A99, L149, and A159 formed hydrophobic interactions. We proposed that dasatinib might be modified at R1, R2, and R3 to improve its inhibitory efficacy against SIK2. Replacement of R1 with a six-membered ring, R2 with an aromatic ring, and R3 with a polar group can enhance the ligand binding affinity for SIK2.

Clarification of dasatinib–SIK2 binding provided insights into the inhibitory mechanism of this drug at the atomic level. This discovery may help guide the design of novel SIK2 inhibitors. However, the present study did not empirically validate the computed and modeled results. Dasatinib–SIK2 binding offered a perspective into inhibitor/SIK2 binding but the other inhibitors may bind dasatinib via different mechanisms. In future studies, we will investigate the mechanism by which other inhibitors such as HG-9-91-01, KIN112, MRT67307, and MRT199665 bind with SIK2. The T-loop conformation may be indicative of the regulation of the activation of SIK2. Further, we will examine various states of the T-loop conformation to elucidate the mechanism by which it influences inhibitor–SIK2 binding.