Molecular dynamics study on inhibition mechanism of CDK-2 and GSK-3β by CHEMBL272026 molecule

CDK-2 and GSK-3β are kinases that exhibit a high level of similarity in the context of amino acid sequence and structural properties of the active site. Both considered enzymes fulfill an important role in the regulation of cell cycle, and disorders in their proper functioning can trigger many serious diseases, even cancer. Identification of potential inhibitors with selective properties is not a trivial task. The potential group of compounds with such properties comprises indirubin derivatives and their analogs. The ligand molecule was docked to two enzymes, namely CDK-2 and GSK-3β, with the use of AutoDock Vina package. The stability of obtained complexes was evaluated with the use of molecular dynamics simulations realized with the use of AMBER11. Both complexes obtained during docking stage are stabilized by network of hydrogen bonds containing similar number of interactions. However, the time evolution of these systems shows significant discrepancies in terms of stability of these interactions and conformational flexibility of ligand molecule. Differences observed for dynamics and structural properties also confirm the values of binding affinity for examined systems. CHEMBL272026 molecule exhibits binding properties toward both considered kinases. However, binding affinity values and examined interactions denoted for both complexes show significant discrepancies in stability of both systems, indicating the higher inhibiting properties toward CDK-2 active site. The discovery of inhibitors with selective potential is an important task, and investigating the mechanism of kinases inhibition may reveal factors responsible for the increase in selective properties what can contribute to the development of new groups of compounds.


Introduction
Kinases are a large group of enzymes responsible for regulation and control of most cellular processes. Many of them, even though they are involved in regulation of different cellular processes, exhibit substantial similarity in terms of the structure of ATP binding pocket. Such a situation may substantially impede the discovery of selective drugs dedicated as inhibitors toward specific biological target. Cyclin-dependent kinase 2 (CDK-2) [1] and glycogen synthase kinase-3 (GSK-3b) [2] are an example of enzymes with very well-known structure of active site. Close resemblance of their amino acid sequence is quite significant (25 % identity and 41 % similarity) [3,4] what contributes to many problems with selective inhibition. The first of considered kinases, namely CDK-2, is serine/threonine kinase which plays a crucial role in regulation of the proliferation of eukaryotic cells. The crucial role in the control of the cell division cycle is the reason why CDKs are one of the most important therapeutic targets in anticancer drug research [5][6][7]. The most important activities of cyclin-dependent kinase 2 (CDK-2) are observed during the G1 phase of the cell cycle, and disturbances in proper activity of this enzyme can be associated with many diseases, even with cancer [8][9][10][11][12].
One of the most important aspects of drug discovery is searching for compounds with high potential of selectivity toward examined biological targets. Recent investigations show a large group of ATP-competitive inhibitors of CDK-2 and GSK-3b which exhibit similar affinity toward both considered kinases [21,22]; however, there are also known groups of compounds with identified selectivity potential toward examined kinases like indirubin derivatives and analogs [23][24][25][26][27][28][29][30]. Investigations on the mechanism of kinases inhibition may reveal factors responsible for increase in selective properties and contribute to the development of new groups of compounds.

Docking procedure
The structure of ligand molecule CHEMBL272026 ((3Z)-N-[3-(1H-Imidazol-1-yl)propyl]-2-oxo-3-[(4-sulfamoylphenyl) hydrazono]-5-indolinecarboxamide) was downloaded from CHEMBL database, while both proteins used during docking procedure, namely CDK-2 (PDB ID 1E9H) [1] and GSK-3b (PDB ID: 1Q41) [2], were obtained from Brookhaven Protein Database PDB. Docking stage of computations was realized with the use of AutoDock Vina [31] utilizing united-atom scoring function. During initial preparations, all water molecules were removed from crystal structures of both enzymes. All structures used during docking, namely ligand and proteins, contained only polar hydrogen atoms. All preparation steps were realized using AutoDock Tools package. The search space for both active sites was limited by the following dimensions: 16 9 14 9 22 for CDK-2 and 22 9 22 9 16 for GSK-3b, and all calculations were realized with exhaustiveness parameter equal to 25. Docking procedure, in the case of both active sites, was repeated ten times. To the next stage of simulations were selected the ligand conformations, which provide occurrence of key interactions, characteristic for complexes of indirubin derivatives with active sites of CDK-2 and GSK-3b.

Molecular dynamics simulation procedure
During molecular dynamic simulations, the time evolutions of two complexes created by CHEMBL272026 with CDK-2 and GSK-3b protein subunits were analyzed. The ligand structure was characterized using the AMBER force field parameters, and the atomic charges were calculated according to the Merz-Kollmann scheme via the RESP procedure [32] at HF/6-31G* level of theory. Each system was neutralized and immersed in a periodic TIP3P water box (GSK-3b: 97 9 90 9 90 Å , 19,953 water molecules; CDK-2: 85 9 92 9 76 Å , 14,817 water molecules). Both systems after two stages of minimization were heated up to 300 K by 100 ps of initial MD simulation. The temperature of the system was controlled by the Langevin thermostat [33]. The periodic boundary conditions and SHAKE algorithm [34] were applied to 70 ns of molecular dynamic simulation. For both considered complexes, the first 10 ns of the simulation time was treated as equilibration stage, and the next 60 ns of trajectory was used in the analysis of dynamic and structural properties of considered complexes. The evaluation of binding affinity of the ligand toward the active site was conducted with the use of the molecular mechanic/Poisson-Boltzmann surface area (MMPBSA) method [35]. In all molecular dynamics simulations, AMBER 11 package [36] was used with ff99SB Force Field [37]. All structural analyses including hydrogen bonds, distances, dihedral angles and root-mean-square deviation (RMSD) of atomic position were performed by the VMD package [38]. For identification of hydrogen bonds were used the following criteria: distance between donor (D) and acceptor (A) \ 3.5 Å , distance H-A \ 3 Å and angle D-H-A [ 90°. All RMSD calculations performed for ligand molecule and enzymes were performed relative to the initial structure of complexes obtained during the docking stage; in these calculations, the hydrogen atoms were excluded.

Results
Structural and energetic properties of complexes obtained during docking stage Both active sites of enzymes examined in this paper exhibit a high level of similarity. Previous studies allowed for identification of the most important amino acids which should be involved in the creation of bindings with potential inhibitors [1,2,39]. Their localization in the conformational space of binding pocket is quite similar for both proteins. In the case of glycogen synthase kinases, the most important role is fulfilled by aspartic acid (ASP 133) and valine (VAL 135). These subunits are involved in creation of one and two hydrogen bonds [2,39], respectively, while for the cyclin-dependent kinases, such role is played by guanine (GLU 81) and leucine (LEU 83): The first one creates one hydrogen bond, while the second one provides two bindings [1,2].
The docking simulations confirm that also in the case of interactions with CHEMBL272026 molecule, these amino acids play a significant role in stabilization of the complexes. The structure of the ligand molecule is presented in Fig. 1, marked atoms are involved in creation of bindings with protein active site, as identified during docking procedure. Complex with CDK-2 is stabilized by 6 hydrogen bonds, namely Ligand Mutual orientation of hydrogen atom H3 from ligand and carboxyl oxygen from LEU 83 allows for interaction; however, the measured distance exceeds the value adopted as the limit for the effective hydrogen bond interaction. In the case of the complex of the CHEMBL272026 ligand with GSK-3b, 7 hydrogen bonds were found: Ligand The average binding affinity of the considered ligand toward both active sites adopts very similar values, in the case of CDK-2 -9.06 (±0.13) kcal/mol and -8.68 (±0.27) kcal/mol for GSK-3b. Values recorded for both complexes indicate a slightly higher affinity of CHEMBL272026 toward CDK-2 protein.

Stability of complexes during molecular dynamics simulation, RMSD analysis
Structural and dynamic properties of the most representative structures obtained during docking procedure were evaluated with the use of molecular dynamics simulations. Structural stability of complexes, during whole simulations, was described with the use of RMSD values for ligands and proteins.
Averaged values presented in Table 1 and graphs attached in Fig. 3 unambiguously indicate the high stability of analyzed systems. In the case of all distributions describing ligands and both proteins, the first 10 ns of trajectory is a sufficient equilibration period; in the case of the rest of conformations, we can observe fully equilibrated structures without any important fluctuations, indicating the dynamic instability of systems. The differences in averaged values and standard deviations characterizing ligands in both systems and also distributions observed on the graphs (Fig. 3a, c) suggest a higher flexibility and mobility of the ligand in the complex with GSK-3b active site compared with the second analyzed protein complex. In the case of considered proteins, also some differences are observed alike for averaged values, standard deviations and distributions on graphs. All of them clearly indicate a higher structural flexibility of GSK-3b than CDK-2 protein.
Evolution of conformational properties of CHEMBL272026 structure during molecular dynamics simulations During molecular dynamics simulations, both structures of complexes evolved and structural changes in terms of conformational properties of ligands and hydrogen bonds were observed. Structural changes of ligand molecules are the result of flexibility of both side chains connected to indoline core of the molecule. The conformational changes in ligand structure in complex with CDK-2 active site are defined with the use of three dihedral angles formed by the following atoms: C2-C10-C20-N7; C21-C18-N1-N4; and C3-C14-S1-N3. All angles are defined in Fig. 4a. The dihedral angle describing orientation of amid group relative to indoline ring (C2-C10-C20-N7) in the initial structure takes a value -35.14 deg, but during molecular dynamic simulation, it is shifted to the positive range of values what can be observed in Fig. 4e. Such conformational change could be caused by two factors. One of them is the creation of internal hydrogen bond between peptide hydrogen atom H1 with nitrogen N2 from imidazole ring. Such interaction is observed during the whole simulation time. This coincidence is also confirmed by the graph presented in Fig. 5a.
Mutual relation between hydrogen bond length and values of considered dihedral angle clearly shows that increase in angle value contributes to the shortening of hydrogen bond length. The lengths of all considered hydrogen bonds indicate the weak character of interactions; however, it may be an important contribution in the context of  Fig. 4f. Such structural change is probably determined by the creation of internal hydrogen bond between hydrogen H3 and nitrogen N2 from imidazole ring. The mutual relation between occurring hydrogen bond and changes in values of C21-C18-N1-N4 dihedral angle is presented in Fig. 5c. In the presented population of hydrogen bonds, weak interactions are dominant.
Next conformational change is related to orientation of sulfamide group relative to aromatic ring defined by C3-C14-S1-N3 dihedral angle. Among the structures collected during the molecular dynamics simulation, dominate con Conformational fluctuations of CHEMBL272026 molecule in complex with GSK-3b show distributions of values of two dihedral angles defined by the following atoms: N1-N4-C12-C7 and N3-S1-C14-C4. In the case of the first presented angle, initial value (112) characterizing structure obtained during docking stage shifts to higher values equal nearly to 180 deg. Such conformational modification changes the orientation of aromatic ring relative to imidazole core of ligand. In the case of sulfonamide group, similar modification is observed like the one denoted for structure complexed with CDK-2 protein. The N3-S1-C14-C4 angle value decreased from -166.14 deg to values placed in the scope from -60 to -110. This conformational shift destabilizes two hydrogen bonds identified during docking stage. In this case, any new hydrogen bond interaction which could stabilize this part of ligand was not found.   Table 2 show that this binding can be classified as medium strength interaction (90 % of population), and its length fluctuates mainly in the range from 1.7 to 2 Å . Second binding, created by atoms from indoline core of ligand, is localized between oxygen atom O3 and peptide hydrogen atom of leucine (Ligand (O3) … (HN) LEU 83). This binding is also observed during all simulation time. However, the share of medium strength interactions is reduced to 65 %, and the rest of conformations create weak hydrogen bonds. Next important hydrogen bond is localized in the closest vicinity of indoline part of ligand. Here, oxygen atom O4 from amide group creates bindings with hydrogens from amino group of lysine (LYS 33). The rotation of amino group during molecular dynamics simulation caused all three hydrogens to participate in the creation of bindings; therefore, for presentation of these interactions, values of distance between oxygen O4 atom and nitrogen atom from lysine amino group (Ligand (O4) … (N) LYS 33) were used. Values presented in Table 2 show that these interactions are observed almost for all conformations collected during molecular dynamic simulation; taking into account individual data for each interaction, it can be stated that the biggest fraction represents medium strength hydrogen bonds (70 %). The atoms of sulfamide group also participate in stabilization of complex with CDK-2 active site. Hydrogen bonds identified during docking stage disappear and are replaced by new ones created with aspartic acid (ASP 86). First interaction is created by peptide hydrogen atom of ASP 86 with nitrogen N3 from sulfamide group (Ligand (N3) … (H) ASP 86). These bindings are found in 92 % of conformations collected during molecular dynamics simulation, but only less than 20 % of them are medium strength interactions, while the rest of them are weak hydrogen bonds. Last important hydrogen bond observed between CHEMBL272026 ligand and CDK-2 protein is created by hydrogen atom H21 from sulfamide group and oxygens O1 and O2 atoms from carboxyl group of aspartic acid (ASP 86). Rotation of carboxyl group caused that both oxygens create interactions with hydrogen atom H21, but the dominant role in these interactions is played by oxygen atom O1. The presence of such hydrogen bonds was found in 92 % of conformations. Values presented in Table 2 clearly indicate that they are mainly interactions of medium strength with weak hydrogen bonds represented by less than 30 % of population.
In the active site of CDK-2 enzyme, two phenylalanines, namely PHE 80 and PHE 82, are localized. Previous studies indicate that their aromatic rings can play significant roles in the creation of hydrophobic interactions [31]. Figure 6 shows values describing distribution of distances between aromatic rings of phenylalanines, aromatic ring of indoline and phenyl ring from side chain of ligand. In the b Fig. 4 Definitions and values of dihedral angles describing structural properties of CHEMBL272026 ligand molecule obtained during docking procedure (a, c) and representative structures from molecular dynamic simulation (b, d). Distributions presented on graphs (ei) describe fluctuations of values for individual dihedral angles defined for ligand molecule during molecular dynamics simulation  case of phenylalanine PHE 80, the distance and orientation of aromatic system relative to indoline for more than 90 % of conformations collected during molecular dynamics simulations allow for the presence of stacking interactions between aromatic systems. The second presented distribution, describing mutual orientation of phenylalanine PHE 82 and phenyl ring, unambiguously indicate the lower affinity of considered systems. The higher distances between aromatic systems caused that only in the case of 50 % of conformations stacking interaction could be considered.

Interactions stabilizing GSK-3b complex
The complex of CHEMBL272026 ligand with GSK-3b protein is stabilized by four hydrogen bonds. The indoline core of ligand creates two hydrogen bonds. First of them is localized between hydrogen H2 atom and carbonyl oxygen of aspartic acid ASP 133 (Ligand (H2) … (O) ASP 133). Data presented in Table 2

Binding affinity analysis
The energetic properties of CHEMBL272026 complexes with CDK-2 and GSK-3b proteins were characterized with the use of the molecular mechanic/Poisson-Boltzmann In the case of the structure resulting from the complex with CDK-2, structural changes promoting creation of internal hydrogen bonds were observed which stabilized ligand structure and restrict the side chains mobility of considered molecule. The observed structural stiffness of ligand from complexes with CDK-2 protein may be one of the factors increasing its binding affinity toward protein active site. The analysis of molecular dynamics trajectories obtained for both complexes clearly show significant differences in the stability of hydrogen bonds of examined complexes. In the case of both complexes, disappearance of part of interactions identified during docking stage were observed. However, the changes observed for complex with CDK-2 are related with replacing of some interactions by another more stable bindings (Ligand (O2) … (HN) GLN 85, Ligand (HN3) … (O) HIE 84, Ligand (O2) … (HN) LYS 89 replaced by Ligand (N3) … (H) ASP 86, Ligand (H21)… (O1) ASP 86 and Ligand (H21)… (O2) ASP 86), while in the case of complex with GSK-3b, interactions created by sulfamide group disappear completely. In the case of both complexes, interactions created by indoline core and amide oxygen atom O4 have similar characteristics. During studying of molecular dynamics, an additional aspect was found that influences the CDK-2 complex stability, namely the possibility of appearance of stacking interactions between phenylalanine and aromatic systems of ligand. Such interactions could increase binding affinity of ligand toward active site of enzyme.
The observed differences in the context of structural stability of ligand molecule, quantity and diversity of interactions correlate with the global values of Gibbs free energy evaluated for considered complexes. Twofold lower affinity of ligand toward protein reported for the complex with GSK-3b indicates significant differences in inhibition properties of CHEMBL272026 molecule. Such important discrepancies reported in binding process may indicate potential selectivity of considered ligand in interactions with CDK-2 and GSK-3 b proteins.