Molecular dynamics simulation studies of GSK-3β ATP competitive inhibitors: understanding the factors contributing to selectivity

Glycogen synthase kinase-3 is a constitutively acting, multifunctional serine threonine kinase, the role of which has been implicated in several physiological pathways and has emerged as a promising target for the treatment of type-II diabetes and Alzheimer’s disease. In order to provide a detailed understanding of the origin of selectivity determinants of ATP competitive inhibitors, molecular dynamics simulations in combination with MM-PBSA binding energy calculations were performed using crystal structures of GSK-3β and CDK-2 in complex with 12 ATP competitive inhibitors. Analysis of energy contributions indicate that electrostatic interaction energy dictates the selectivity of ATP competitive inhibitors against CDK-2. Key interactions as well as residues that potentially make a major contribution to the binding free energy were identified at the ATP binding site. This analysis stresses the need for the inhibitors to interact with Lys85, Thr138, and Arg141 in the binding site of GSK-3β to show selectivity. The residue-wise energy decomposition analysis further suggested the additional role of Gln185 in determining the selectivity of maleimides. The results obtained in this study can be utilized to design new selective GSK-3 ATP competitive inhibitors.

The advancement of computational power and computational methods like molecular docking, 3D-QSAR, pharmacophore mapping, virtual screening, and molecular dynamics (MD) simulations have provided breakthroughs in the field of drug discovery. Molecular docking methodology helps to predict good binding poses for the ligands but generally fails to predict binding affinities in good correlation, probably due to the evasion of protein flexibility (Ferrara, Gohlke, Price, Klebe, & Brooks, 2004;Srivani, Srinivas, Raghu, & Sastry, 2007;Warren et al., 2006). The generation of QSAR and pharmacophoric model is based on the features of known ligands and hence can effectively be applied to design new ligands which belong to same or similar class (Akamatsu, 2002;Kubinyi, Folkers, Martin, & Martin, 2002). The reason for this discrepancy is negligence of important features of active site and protein flexibility, which varies from inhibitor to inhibitor. The development of methods like MD simulation Monte Carlo simulations, free energy perturbation (FEP), (Zwanzig, 1954) thermodynamic integration (TI), (Kirkwood, 1935) linear interaction energy analysis, (Åqvist, Medina, & Samuelsson, 1994), and MM-PBSA/ GBSA (Kollman et al., 2000;Srinivasan, Miller, Kollman, & Case, 1998) calculations provide very efficient key to the above-mentioned problems. TI and FEP are more rigorous and accurate methods to calculate the binding free energy but are computationally expensive and time consuming in contrast to MM-PBSA and MM-GBSA methods which are faster by several magnitudes (Ajay and Murcko, 1995;Gilson, Given, Bush, & McCammon, 1997;Jorgensen, 1989;Kollman, 1993;Simonson, Archontis, & Karplus, 2002).
Considering the importance of GSK-3 which is a well-established target for drug discovery, MD simulations have been previously performed to address the issues like (i) the preferential binding of enzyme to the phosphorylated substrate than non-phosphorylated substrates, unique P + 4 primed phosphorylation specificity of GSK-3β for its substrates (Lu, Zou, et al., 2011), (ii) the role of conserved water molecules in the crystal structures of GSK-3β (Hayes et al., 2011;Lu, Lv, et al., 2011), (iii) allosteric regulation of GSK-3β (Buch et al., 2011), (iv) the loss of catalytic activity due to mutation (Sun, Jiang, Yu, Luo, & Zou, 2008;Zhang, Jiang, Zou, Yu, & Zhao, 2009), (v) peptide and GSK-3 binding mode analysis (Tang et al., 2011), (vi) the inhibition of GSK-3 by metals (Lu et al., 2013;Sun, Jiang, Yu, Luo, & Zou, 2011), (vii) binding mode analysis of maleimide with GSK-3 (Jena, 2012), (viii) TI MD simulation study to evaluate binding free energy associated with GSK-3β (Lee, Hsu, Liu, Hsu, & Sun, 2014), (ix) auto inhibition of GSK-3β by Arg4 and Arg96 (Mou et al., 2014), etc. MD simulations have also been carried out to address the issue of selective inhibition of GSK-3. Pande and Ramos explored the structural basis for the GSK-3β binding affinity and selectivity against CDK-2 (Pande & Ramos, 2005). Chen et al. studied the selective inhibition of GSK-3β by paullones (Chen, Cui, Cheng, Zhang, & Ji, 2011). Zhang et al. evaluated the structural features responsible for selective inhibition of GSK-3 by dibromocantharelline (Zhang, Zhong, Yan, & Jiang, 2011). Rastelli and co-workers assessed the protein kinase selectivity of several kinases including GSK-3 (Muzzioli, Del Rio, & Rastelli, 2011). However, the MD simulation studies related to the selectivity of ATP competitive inhibitors have been mainly limited to a specific class and moreover, these reported simulation studies have not taken into account the thermodynamic viewpoint and have not provided the detailed understanding of binding determinants. Therefore, in this work, an exhaustive and systematic MD simulation study has been taken up to address the selectivity issue of ATP competitive inhibitors from the thermodynamic perspective. Complexes of 12 inhibitors with GSK-3β and CDK-2 belonging to 5 different classes were selected for MD simulation and energy component analysis studies to obtain the information on physiochemical forces that drives the selective inhibition GSK-3β. The binding mode and interaction profile of the inhibitors against GSK-3β as well as CDK-2 have also been explored and compared to pinpoint the selectivity determinants. On the basis of obtained results, we present the systematic discussion on the energy components and structural determinants.

Molecular docking
The molecular docking was carried out to obtain the starting structures of the inhibitor-protein complexes for which the co-crystallized structures were not available.
The three-dimensional (3D) structures of ligands were built using sketch module of SYBYL7.1 and subjected to 1000 cycle minimization with standard Tripos force field and .005 kcal/mol energy gradient convergence criterion using Powell's method (Powell, 1964). 3D structures of these compounds were then prepared using LigPrep module of maestro implementing OPLS_2005 force field and ionic states (carboxylates) for the ligands at pH values of 7.0 ± 2.0. Molecular docking was performed using Glide software (Grid-based Ligand Docking with Energetics), (Glide, 2012) with the standard precision mode. Initially, to check the reliability of docking protocol, the bound ligands (as given in X-ray crystal structures) were redocked (see Supplementary Information) in the active sites of respective enzymes and RMSD values of the docked poses were compared with the X-ray crystallographic conformation of the ligand. Before performing the docking experiments, conformational search was performed using confgen module of maestro. Top five conformations of low energy were docked into the ATP binding site of the GSK-3β and CDK-2.

MD simulation
AMBER 11.0 program was used for MD simulations of selected starting structures (Case et al., 2005). The force field "leaprc.gaff" (generalized amber force field) was used to prepare the ligands, while "leaprc.ff03" was used for the enzyme. Each system was placed in an octahedral box (with a 10.0-Å boundary) of TIP3P water and chloride ions were added to neutralize the system. Each complex was subjected to minimization (500 steps of steepest descent followed by 500 steps of conjugate gradient method), 50 ps of heating, and 50 ps of density run with weak restraints on the complex followed by 4-ns constant pressure equilibration at 300 K. A cut-off of 8.0 Å was used for MD simulations. All long-range electrostatics were included by means of a particle mesh Ewald method. All hydrogen-heavy atom bonds were constrained by the SHAKE method, and simulations were performed with a 2-fs time step. Langevin thermostat was used to maintain the temperature of the system at 300 K. The same conditions as the final phase of equilibration were used for production run, and the coordinates were recorded every 10 ps. The periodic boundary conditions were used during MD simulations. Before submitting for the production run of 24 ns, root mean squared deviation plots (RMSD) were constructed to make sure that system has equilibrated.

MM-PBSA calculations
One thousand and four hundred equally spaced snapshots of each complex (every 10 ps) were generated from the MD trajectories of last 14 ns of MD simulation run. All water molecules and counter ions were removed before binding energy and residue-wise energy decomposition calculations, which was done using python script of AMBER 11.0 program.
The DG Bind-PB between ligand and receptor to form a complex can be calculated as follows.
where ΔE MM is the total gas-phase energy, DG SOLðPBÞ is the sum of non-polar and polar contributions to solvation calculated by PB. DE internal is the internal energy arising from bond angle and dihedral terms in the MM force field and is zero in case of single trajectory approach. DE electrostatic is electrostatic energy as calculated by the molecular mechanics force field. DE vdw is van der Waals contribution from MM. ΔG PB is the non-polar contribution to the solvent free energy and was estimated by SASA determined using a water probe radius of 1.4 Å. The surface tension constant γ was set to be .00542 kcal/mol/Å 2 . DG SAðPBÞ is the electrostatic contribution to the solvation energy calculated by the PB method. Dielectric constants for solute and solvent were set to be 1 and 80, respectively. Vibrational entropy contributions can be estimated by classical statistical thermodynamics, using normal mode analysis. As our goal is not to obtain the absolute Gibbs free energy but to analyze the comprehensive interaction features, the entropy contribution was not included in the study due to high computational cost.

Inhibitor-residue interaction decomposition
The interactions between each residue and inhibitor were computed using the MM-PBSA decomposition process applied in the MM-PBSA module of AMBER11.0. The binding interaction of each inhibitor-residue pair includes three terms, van der Waals contribution (ΔG vdW ), electrostatic contribution, and solvent contribution (ΔG SA(PB) ).
where DG vdw and DG elc are non-bonded van der Waals interactions and electrostatic interactions between the inhibitor and each GSK-3β residue, which can be computed using the python script of AMBER11.0. The ΔG PB of solvation was computed using Poisson-Boltzmann model.

Structural analysis and multiple sequence alignment
The crystal structure analysis showed that Asp133 and Val135 form hydrogen bond interaction with inhibitor and is observed in almost all the co-crystals of ATP competitive inhibitors with GSK-3β. The other polar interaction observed in the co-crystals of GSK-3β with ATP competitive inhibitors is water mediated which is between Thr138, Gln185, and inhibitors. In case of 1Q4L, an electrostatic interaction is observed between carboxylic group of ligand and Arg141. The hydrophobic pocket in ATP binding domain of GSK-3β is formed by Ile62, Val70, Ala83 on the top and Leu188 at the bottom. The complementary part of the hinge segment consists of Leu132-Asp133-Tyr134-Val135-Pro136 and forms one portion of the pocket, while Arg141 forms another part of the pocket by orienting the hydrophobic side chain toward it. Therefore, a sequence alignment was performed to identify the difference in the ATP binding pocket of GSK-3β and CDK-2, which can then be analyzed in terms of binding free energy contribution to derive the useful information, which can further be utilized to design potentially selective GSK-3β ATP competitive inhibitors. The sequence alignment (Figure 3(A)) shows that top and the bottom part of the hydrophobic pocket in GSK-3β and CDK-2 are formed of identical amino acids (Ile10, Val18, Ala31, and Leu134 in CDK-2), the hinge segment in CDK-2 is similar in nature (Phe80-Glu81-Phe82-Leu83-Hie84 for CDK-2 and Leu132-Asp133-Tyr134-Val135-Pro136 for GSK-3β) and in place of Arg141, there exists Lys89 which again is similar in nature (Figure 3(B)-(D)).

Molecular docking
Molecular docking simulations were performed to obtain the starting pose of the ligands in complex with GSK-3β and CDK-2 for which the crystal structures were not available. A maximum of 10 docking poses per ligand were generated in each case and analyzed further for the binding mode and intermolecular interactions. The cocrystallized complex for ligands L1, L4, L8, L10, L11, and L12 were available and hence were obtained from protein data bank (PDB code 1Q4L, 4ACC, 1Q5K, 1Q3W, 1Q41, and 1UV5, respectively) to be used for the MD simulation study. The starting pose for other ligand in complexation with GSK-3 was obtained through molecular docking methodology using PDB code 1Q4L (L2 and L3), 4ACC (L5-L7), and 1Q3W (L9). Similarly, the co-crystallized structure for L4 in complexation with CDK-2 was obtained from PDB used for the MD simulation, while for other ligands, the CDK-2 complexes were obtained through molecular  (4ACM)] or through manual modeling (L1, L2, and L8, by aligning the co-crystallized complex of GSK-3β with CDK-2 and merging the ligands into the ATP binding cavity of CDK-2) using PyMOL v0.99.

MD simulation
MD simulations were carried out for 24 ns to ensure that equilibration state of the complex lasted enough for further binding energy analysis. The convergence of all the simulations was analyzed in terms of energy components and RMSD from initial structure.

RMSD examination
To explore the dynamic stability of 24 protein/ligand complexes, root mean square deviation (RMSD) values for the heavy atoms of the GSK-3β and CDK-2, during the production phase, relative to the starting structure and first frame obtained after 4-ns simulation run were calculated (since in some cases, desired starting poses were not obtained with the help of docking methods, and hence were manually modeled). The RMSD plots indicate that out of the 24 protein/ligand complexes submitted for MD simulations, 23 protein/ligand complexes achieved good equilibrium within 2 ns (Figure 4). The trajectory for complex of CDK-2 with ligand L7 showed larger fluctuations and was not able to achieve good equilibrium even after the simulation run of 20 ns (Figure 4(B)).

Energy examination
To calculate reliable binding energy values, the average values calculated by the MM-PB/GBSA must be converged. To assess the convergence of binding energy, cumulative mean (μ) and standard error (σ) of binding energy as a function of forward time interval from the beginning of simulation run were calculated (see Supplementary Information). The mean values, converge to within 1.5 kcal/mol at t = 10 ns of the trajectory in the forward direction for most of the systems, which is consistent with the RMSD analysis. The averaged energy values are converged which is further supported by the low value of standard error (.3) for the mean energies calculated for the last 14 ns. The convergence of binding free energy at t = 10 ns indicates that simulation run of more than 10 ns is required to get the reliable estimates of binding energy and therefore the snapshot generated during the simulation run of last 14 ns have been used for the purpose of analysis. The binding energy calculated from snapshot extracted over the simulation run of last 14 ns was also examined to estimate the binding efficiency of ligands in GSK-3 and CDK-2. The calculated binding energy was also correlated with the reported IC 50 values (as K i values were not reported for all the compounds) and it was found that calculated binding energy and experimentally observed affinity follow the same trend within a class, i.e. binding affinity of L1-GSK-3β = L2-GSK-3β > L3-GSK-3β. Similarly, the trend was consistent with respect to the binding energy of Ligand-GSK-3β vs. Ligand-CDK-2 within a class, i.e. L1-GSK-3β > L1-CDK-2, L2-GSK-3β > L2-CDK-2 and L3-GSK-3β ≈ L3-CDK-2 (see Supplementary Table S7). Careful examination of relative binding energy showed that GSK-3 selective inhibitors form more stable binding with the GSK-3β as compared to the CDK-2, while the non-selective ligands either showed higher binding energy values for the CDK-2 or binding energy values were comparable for both the target enzymes (consistent with the experimentally observed affinity of ligands with respect to the target enzymes). In order to elucidate the binding mechanism of selective and non-selective ligands, individual contributions to the binding energy were also examined and it was observed that the van der Waals interactions between the binding partner, the nonpolar interactions with the solvent including the contribution from the hydrophobic effect and the electrostatic interactions played the dominant role in determining the binding of the selective inhibitors towards GSK-3β. The most stabilizing value for electrostatic interactions was found for ligands L1 and L2 complexed with GSK-3β ( Figure 5, Table 2). Intermolecular electrostatic interactions include contribution from direct hydrogen bonds, water-mediated hydrogen bonds, and other polar interactions between protein and ligand. In case of GSK-3β, the hydrogen bonds are generally formed by the inhibitors with Asp133 and Val135, while in case of ligands L1 and L2, high amount of electrostatic energy is contributed by interaction of carboxylic group of the ligand (L1 and L2) with the positively charged Arg141, Thr138, and Gln185 located at the opening of the cavity. The association of inhibitors to the GSK-3β and CDK-2 is opposed by an unfavorable desolvation of polar groups yielding a contribution of 74-90 kcal/mol. This unfavorable desolvation of polar groups is only partially compensated by favorable electrostatic interactions in most of the enzyme inhibitor complexes except in case of L1 and L2 where it is favored for GSK-3β association by 3 kcal/mol.
Anilinomaleimide (L1, L2) vs. indolylmaleimide (L3) As expected, for the GSK-3β complexed with ligands L1 and L2, the binding energy (40.99 and 38.41 kcal/mol) was found to be more favorable than binding energy of L1 and L2 complexed with CDK-2 (18.87 and 20.87 kcal/mol respectively), while binding energy of L3 complexed with GSK-3β and CDK-2 was found to be comparable in strength (26.06 and 23.81 kcal/mol, respectively). For ligands L1 and L2 complexed with GSK-3β, the gas-phase van der Waals interaction energy is favored by 4.9 and 7.6 kcal/mol as compared to their CDK-2 counterpart. For ligand L3 complexed with GSK-3, the van der Waals interaction energy is favored for GSK-3 by 8.6 kcal/mol. The gas-phase electrostatic interaction energy of L1 and L2 also favors binding with GSK-3β by 109.82 and 116.68 kcal/mol, respectively, while for L3 the electrostatic interaction energy is .41 kcal/mol in favor of CDK-2. The difference in electrostatic energy pattern can be attributed to the presence of carboxylic group in ligands L1 and L2, which form salt bridge with the Arg141, while Lys at the similar position in CDK-2 does not form such strong electrostatic interaction. In case of L3, the indole part of the ligand, which possess a hydroxyl group and orients toward the Arg141, and therefore does not make any significant electrostatic interactions with GSK-3β. Considering only the gas-phase interaction energies for the maleimide series of ligands, it can be clearly stated that the electrostatic interaction energy plays a dominant role in determining the selectivity. Also, when solvent-phase energies are analyzed, it is observed that in case of L1 and L2, the electrostatic interaction energy favors the binding with GSK-3 by 17.13 and 9.85 kcal/mol and for the L3, the binding is opposed by 6.64 kcal/mol. The solvent-phase van der Waals interaction energy for L1, L2 and L3 favors the binding with GSK-3β by 4.97, 7.96, and 8.88 kcal/mol, respectively. From the analysis of gas-phase and solvent-phase energies, it can be stated that the electrostatic energy component plays a dominant role in determining the selectivity for the maleimide class of compounds.
Pyrazines (L6, L5 vs. L4, L7) For pyrazine derivatives, the binding energy of ligands L5, L6 (selective), L4, and L7 (non-selective) for GSK-3 (CDK-2) are 29. 13, 30.24, 28.63, and 24.65 kcal/mol (22.70, 26.00, 25.66, and 24.23 kcal/mol), respectively. The gas-phase energy component analysis shows that electrostatic interaction energy in case of L5 and L6 favors binding with GSK-3β by 9.51 and 13.87 kcal/mol, respectively, against CDK-2, while for L4 and L7 the gas-phase electrostatic interaction energy favors binding with CDK-2 by .01 and 5.00 kcal/mol respectively. The gas-phase van der Waals interaction energy for L5 favors GSK-3β binding by 3.88, while for L6 van der Waals interaction energy favors the binding to CDK-2 by .36 kcal/mol. The van der Waals interaction energy for L4 favors the binding to CDK-2 by 1.35 and for L7 it favors GSK-3 by 5.56 kcal/mol. However, the analysis of solvent-phase energy contributions indicates that electrostatic interaction energy favors the association of L4, L5, and L6 to GSK-3β by 4.59, 2.48, and 4.64 kcal/mol, respectively. For ligand L7, the solvent-phase electrostatic interaction energy favors binding to CDK-2 by 5.27 kcal/mol. The solvent-phase van der Waals interaction energy favors the binding of L5 and L7 to GSK-3β by 3.95 and 5.69 kcal/mol, and for ligands L4 and L6, the solvent-phase van der Waals interaction energy favors binding to CDK-2 by 1.62, .40 kcal/mol, respectively. Therefore, in case of L5, both van der Waals and electrostatic interaction energy plays an important role in determining the selective inhibition of GSK-3β, but for ligand L6, it is only the electrostatic interaction which helps in selective inhibition of GSK-3β.

Nitrothiazole-GSK-3 vs. Nitrothiazole-CDK-2 (L8)
Ligand L8 is one of the most selective inhibitor reported till date, the binding energy for GSK-3β and CDK-2 was found to be 28.8 and 26.7 kcal/mol. The gas-phase electrostatic energy was found to be in favor of CDK-2 by 1.14 kcal/mol (26.96 and 28.10 kcal/mol for GSK-3β and CDK-2, respectively). However, in the solvent-phase Figure 5. The bar diagram represents the difference in the binding energy of the GSK-3β inhibitor and CDK-2 inhibitor complex for 12 ligands (L1-L12). The blue color represents the electrostatic contribution and pink color represents the van der Waals contribution to the total binding energy. electrostatic interaction energy is in favor of GSK-3β by~1 kcal/mol due to non-favorable polar solvent contribution toward the binding of L8 with CDK-2 by 2.23 kcal/mol. The gas-phase (36.52, and 35.58 kcal/mol) and solvent-phase (39.20 and 38.23 kcal/mol) van der Waals interaction energies for GSK-3β and CDK-2 are found to be in favor of GSK-3β by~1 kcal/mol.

Paullones (L9 vs. L10)
1-azakenpaullone (L9) is one of the mildly selective inhibitor reported during the early reports of the GSK-3, as per the literature available L9 is 200-fold selective for GSK-3 than for CDK's. The gas-phase and solvent-phase binding energy for ligand L9 favors binding with GSK-3β by 4.3 and 3.2 kcal/mol against CDK-2, while for ligand L10 the binding to CDK-2 is favored by 3.02 kcal/mol in gas phase; however, in solventphase, the binding to GSK-3β is favored by L10 (3.87 kcal/mol). This difference is observed due to the polar solvation energy which was found to be in the favor of GSK-3 by 6.55 kcal/mol. The gas-phase electrostatic interaction energy is found to be unfavorable for binding with GSK-3β by 2.36 kcal/mol for L9, while this unfavorable binding is found to be of 10.20 kcal/mol in strength for L10. The gas-phase van der Waals interaction energy is found to be 6.73 and 7.18 kcal/mol in favor of GSK-3β for L9 and L10, respectively. The nonpolar solvent contribution favors binding with GSK-3β by .15 and .34 for L9 and L10. The above-mentioned analysis of the gas-phase energies shows that the electrostatic interaction energy which opposes the binding of L10 to GSK-3β by 10.20 kcal/mol may have a possible role in determining the selectivity of L9 against CDK-2.
Indirubins (L11 vs. L10) 6-Substituted indirubins are known to be weekly selective for GSK-3β inhibition against CDK-2; therefore, 6-bromoindirubinmonoxime (L11) and indirubinmonoxime (L12) were chosen for this MD simulation study. 6-Bromoindirubin has been reported to be only 16 times selective; however, they are earliest known inhibitors for GSK-3β. The solvent-phase analysis of MM-PBSA binding energy shows that binding of L11 and L12 are favored for GSK-3β by 2.65 and 3.00 kcal/mol, respectively. The gas-phase MM-PBSA binding free energy analysis shows that the binding of L11 to CDK-2 is favorable by 1.18 kcal/mol; however, binding of L12 to GSK-3 is favored by 2.00 kcal/mol. The electrostatic energy favors the binding to CDK-2 by 4.42 and 1.94 kcal/mol (.62 and 1.11 kcal/mol for solvent phase) for L11 and L12, respectively, while the gas-phase van der Waals interaction energy favors the binding to GSK-3β by 3.23 and 3.94 kcal/mol (3.27 and 4.11 kcal/mol for solvent phase).

Per residue inhibitor enzyme interaction
To explore the key binding residues that play an important role in the selective inhibition of GSK-3, the intermolecular interactions between the inhibitors and the protein residues were quantified using pairwise per residue decomposition energy module in AMBER ( Figure 6). The important amino acid residues found to be interacting with inhibitor-GSK-3 (inhibitor-CDK-2) complexes were Ile62 (Ile10), Val70 (Val18), Ala83 (Ala31), Lys85 (Lys33), Val110 (Val64), Leu132 (Phe80), Asp133 (Glu81), Tyr134 (Phe82), Val135 (Leu83), and Leu188 (Leu134). The strength of the interactions varied from 2.62 to 3.98 kcal/mol. It is worth noting that above-identified amino acid residues are important for the molecular recognition, with residues Ile62, Val70, and Ala83 forming the roof of the ATP binding pocket, while Leu188 forming the bottom of the pocket. Residues Lys85 and Val110 form one side of the pocket, hinge binding residues Leu132, Asp133, Tyr134, and Val135 forms the other side of the pocket. Interestingly Asp133, which is one of the hinge binding residues and has been found to form hydrogen bond in many of the reported crystal structures of ATP competitive inhibitors, was found to show unfavorable interaction in case of ligands L1, L2, L9, and L10. The strength of this unfavorable interaction varied from .77 to 1.18 kcal/mol, while for the similar residue in case of CDK-2 (i.e. Glu81), this unfavorable interaction energy varied from .41 to 1.42 kcal/mol. The amino acid residue Arg141 (Lys89) which is located at the opening of the ATP binding cavity was found to show strong electrostatic interaction with L1 and L2 (7.92 and 8.66 kcal/mol, highly selective molecules). For the ligands L5 and L6, which has been reported to show good selectivity against CDK-2, the Arg141 interacts favorably with the inhibitor by 1.54 and 1.01 kcal/mol. Lys85 was found to show favorable interaction with L1 and L2 (~1.5 and 3 kcal/mol, respectively) complexed with GSK-3β, whereas in case of CDK-2, similar residue (Lys33) shows unfavorable interaction (~2 kcal/mol). With respect to marginally selective and weekly selective ligands (L3-L12) complexed with GSK-3β, this amino acid residue does not show any interaction. However, in case of CDK-2 complexed with L4, L7, L8, L11, and L12 (weakly selective ligand), Lys33 is involved in unfavorable interaction of the order of 1.5 kcal/mol. Lescot et al. (2005) through CoMFA analysis observed similar type of interactions in 3-anilinomaleimides with Arg141 and Lys85; however, their role in selectivity was not highlighted. Moreover, Thr138, in case of GSK-3β complex was found to interact favorably, the strongest interaction observed in the case of L1 (2.23 kcal/mol) and weakest observed in case of L7 (.17 kcal/mol). In CDK-2, the amino acid found at the similar position is Asp86, interacts unfavorably with the GSK-3β inhibitors, the most repulsive interaction was observed in case of L1 (2.23 kcal/mol) and the weakest repulsion was observed in the case of L11 (.74 kcal/mol) except L4 where Asp86 interacts favorably by .24 kcal/mol. Besides, the above-mentioned amino acids Gln185 have also been found to show strong interaction with L1 (1.56 kcal/mol) and L2 (3.02 kcal/mol) in case GSK-3β inhibitor complex. In case of L9 and L10, the residue Gln185 shows repulsive interaction of~1.4 kcal/mol. The same repulsive interaction (~1.00 kcal/mol) was also observed in case of L9 and L10, with the residue Gln131 in case of CDK-2-inhibitor complex. In addition to the Gln185, Gly65 and Ser66 have been found to be showing weak favorable interaction with L1 and L2, while Phe67 was observed to show favorable interaction of 1.44 kcal/mol in strength with L3. The above analysis very clearly shows that Lys85 and Thr138 play a very important role in determining the selectivity. This observation further gains strength from the fact that greater the strength of repulsion of residues Lys33 and Asp86 with the small molecule in the CDK-2 complexes, the greater is the selectivity of the molecule for GSK-3β.

Hydrogen bond analysis
Hydrogen bond is the most significant interactions in molecular recognition process. Hydrogen bond occupancy between selected ligand and key residues were analyzed during the last 14 ns of simulation run. The percentage and number of hydrogen bonds were calculated according to the following three criteria: (i) proton donor acceptor distance ≤3.5 Å, (ii) donor-H-acceptor bond angle ≥120°and (iii) hydrogen bond occupancy should be more than 10%. Hydrogen bond occupancy of >75%, 50-75%, and <50% were used to define strong medium and weak hydrogen interaction, respectively. Residue having hydrogen bond occupancy of >100% means that it can form more than one hydrogen bond. The number of hydrogen bonds in case of L1 and L2 with GSK-3β complex was found to be 6 and 5, respectively, whereas a similar ligand in CDK-2 complex formed 3 and 2 hydrogen bond, respectively. Ligand L3 in GSK-3β as well as CDK-2 complexes showed three hydrogen bonds. L4 and L7 in GSK-3β complexes showed two hydrogen bonds, while ligands L5 and L6 the number of hydrogen bond found was three. Ligands L4-L7 complexed with CDK-2 showed three hydrogen bonds. Ligands L8-L10 showed only one hydrogen bond, whereas ligands L11 and L12 showed two hydrogen bonds in GSK-3β as well as CDK-2 complexes. Further analysis showed that ligands L1-L12 are able to form weak to strong hydrogen bonds with Val135 and Leu83 in GSK-3β and CDK-2, respectively (percentage occupancy ranging from 29.59 in case of L3 to 114.58 for L10,  (Lys89) showed percentage occupancy of 11.08 in case of L1. Ligands L5 and L6 were also found to make weak hydrogen bond interaction with Arg141 with percentage occupancy of 12.58 and 16.44, respectively. This is a very interesting observation considering the fact that L1, L2, L5, and L6 are strongly selective inhibitors for GSK-3β. Besides, the above-mentioned hydrogen bonds, ligands L4 and L7 in case of CDK-2 complexes showed weak hydrogen bonds with Lys33 (percentage occupancy of 18.94 and 33.67, respectively). The comparative analysis of energy decomposition across the different classes of compounds suggests that the drug design should focus on increasing the electrostatic interactions with simultaneous increase in hydrophobic interaction. In case of weekly selective ligands, the selectivity is mainly driven by van der Waals energy. The energy decomposition analysis further shows that selectivity in case of L1 and L2 is assisted by decreased polar solvation energy in GSK-3β as compared to CDK-2. For ligands L5 and L6, the electrostatic energy was again found to favor binding with GSK-3β, whereas in case of L4 and L7, the electrostatic energy was found to favor CDK-2, but the polar solvation was found to favor binding with GSK-3β resulting into the very weak selectivity. In case of other (L7-L12) ligands, the selective inhibition is determined by the van der Waals energy and low electrostatic energy. Unfortunately, optimizing the net electrostatic interaction will also lead to the increase in polar solvation energy. Residue-wise energy decomposition analysis indicates that increase in electrostatic energy can be achieved by increasing the interactions with polar residues like Lys85, Arg141, Thr138, and Gln185. Repulsive interaction of Lys33 in case of CDK-2 can be attributed to large size of the ATP binding cavity of GSK-3 as compared to CDK-2. Hydrogen bond analysis showed strong hydrogen bond contacts with the residues Arg141, Thr138, and Gln185 in case of L1 and L2, thus confirming that interaction with the above-mentioned residues will lead to the selective inhibition of GSK-3β. In addition, hydrogen bond contact with Val135 and Leu83 was found to be conserved in all the complexes under study, whereas hydrogen bond with Asp133 was found to be missing in the case of L19 and L10, while it was found to be weak in case of L9.

Conclusions
Twelve inhibitors selected on the basis of selectivity against CDK-2, belonging to five different chemical classes were considered for MD simulation study of 24 ns each. The starting pose for the MD simulation was obtained from protein data bank and molecular docking. Following MD simulations, the RMSD comparison of selective inhibitors in complexes with CDK-2 showed higher RMSD and SD values in comparison to GSK-3inhibitor complexes. MMPBSA calculations were carried out to estimate binding energy and energy component analysis. The energy decomposition analysis of binding energy under gas phase as well as under solvated condition showed that electrostatic interaction energy plays a decisive role in determining the selectivity against CDK-2, as was the case with ligands L1, L2, L5, and L6 (in agreement with the study of Chen et al. reported for selective inhibition of GSK-3 by paullones as against the CDK-5). Residue-wise energy decomposition analysis showed that only 13 amino acids (Ile62, Val70, Ala83, Lys85, Val110, Leu132, Asp133, Thr138, Tyr134, Val135, Arg141, Gln185, and Leu188) for GSK-3, while 9 amino acids (Ile10, Val18, Ala31, Val64, Phe80, Glu81, Phe82, Leu83, and Leu134) for CDK-2 take active participation in enzyme inhibitor interaction. The residue-wise energy decomposition analysis further showed that residue Lys85 (in agreement with the reports of Lescot et al. (2005) for maleimides) and Thr138 interact favorably, while Lys33 and Asp86 (identical residue at similar position in CDK-2) interact unfavorably with the selective inhibitors. Arg141 was also found to interact favorably and strongly in case of highly selective ligand, while Lys89 in CDK-2 was not found to interact with comparable strength in case of selective GSK-3 inhibitors. The residue-wise energy decomposition analysis also showed that Gln185 plays supportive role in selective inhibition of GSK-3 in case of maleimides. Besides, Asp133 was found to show repulsive interaction in the complexes of L8, L9, and L10 with GSK-3β. Therefore, the results of MD simulations clearly points out the role of residues Lys85, Thr138, and Arg141 and electrostatic energy in selective inhibition of GSK-3. The above observations are confirmed through hydrogen bond analysis, which showed L1 and L2 (strongly selective ligands) are involved in strong hydrogen bond interaction with Thr138, Arg141, and Gln185 (% occupancy ≤50). Moreover, in case of L8, L9, and L10 complexed with GSK-3β and CDK-2, Asp133 and Glu81 did not show any hydrogen bond contact (consistent with the results of residue-wise energy decomposition analysis). The present work offers a quantitative and mechanistic rationalization of selectivity for five different structural classes of GSK-3β inhibitors. The result provides the comprehensive analysis of the structure-affinity relationship. Therefore, we anticipate that the current study can give useful insights into the nature of selective inhibitors and may support the future design of more selective GSK-3β inhibitors.