Docking of Glycokinase with Oxo, Sulfo, and Seleno Derivatives of the Carboxamide Activator S41

Inactivation of Glucokinase (GK) is associated with diabetes. Therefore, design of drugs targeting the GK activator site is currently integrated in the strat-egy of the diabetes treatment. The present work investigated the affinity of 30 ligands to GK based on molecular docking using the Gold 5.6 program. Glu-cokinase’s structure was derived from the Protein Data Bank (PDB Code 3S41), while the ligands were seleno, sulfo and oxo derivatives of the co-crystallized carboxamide activator (PDB code: S41). The results of the ligand-protein docking revealed that GK formed thermodynamically stable complexes with all ligands. The main forces stabilizing the complexes are lipophilic interactions, enhanced by hydrogen bonds. Ligand molecular areas responsible for lipophilic and hydrogen bonding contacts with amino acid residues in the allosteric site of GK were evidenced by molecular electrostatic potentials (MEPs). Interestingly, twelve of the S41 derivatives interacted with GK more strongly than the co-crystallized activator, while maintaining the lipophilic contacts with key amino acid residues like Arg63, which are catalytically crucial for therapeutic properties of GK activators (GKAs). It is noteworthy that divalent Se and S atoms were also involved in chalcogen bonds in the GKA site. Those bonds were nearly linear like hydrogen bonds. Such bond directionality should guide the design of pharmacophoric ligands containing chalcogen atoms.


Introduction
Glucokinase (GK) belongs as hexokinase to essential proteins which play a significant role in the animal world [1]. Indeed, glucokinase is one of four hexokinases present in high proportion in some organs such as the liver and pancreas, where it is involved in the hepatic metabolism of glucose and the pancreatic secretion of insulin [2] [3] [4]. The high level of glucose concentration in blood and urines associated with diabetes caused in part by inactivation of the gene encoding glucokinase [2].
Due to weak selectivity and side effects of the commonly used antidiabetic drugs, there is a growing interest in new pharmacological strategies requiring less time and capital to improve the antidiabetic therapy [2]. Disadvantages of current strategies based on experimental screening of molecules to identify any that elicit desired biological response are low rate and assays needing extensive development with validation before use. One of the new approaches is the computer-aided drug discovery (CADD) based on molecular docking in conjunction with molecular dynamic simulation and quantitative structure activity relationship (QSAR). Ligand-protein docking has been used to probe pharmacophoric potential of ligands interacting at the activator site of glucokinase [5]- [14]. A comparative investigation on activator potential of numerous pharmacophores revealed that amide ligands consisting of three flexible cyclic arms joined in Y letter shape showed the highest binding affinity at the GK allosteric site [6].
In the present work, we used the derivatives of N,N-dimethyl-5-(2-methyl-6-((5-methylpyrazin-2-yl)-carbamoyl)benzofuran-4-yloxy)pyrimidine-2-carboxam ide (S41), which show a shape resembling to the Y letter, in order to explore their ability to interact at the allosteric site of GK using molecular docking. The docking was carried out with the GOLD 5.6 program (Genetic Optimization for Ligand Docking) implemented in the Cambridge Structural Database (CSD). Gold fitness scores as well interaction energies were calculated to determine the affinity and the thermodynamic stability of the binding, respectively [15] [16]. The binding affinity of the S41 activator was compared with those of oxo, sulfo and seleno derivatives. GK was derived from the Protein Data Bank (PDB code: 3S41). Das et al. [17] showed that the affinity of a ligand to the GK activa- reported on changes of ligand affinity to N-terminal domain of Heat shock protein 90 (HSP90) when oxygen atoms were replaced by sulfur or selenium atoms [18].

Preparation of 3S41
The crystal structures of the glucokinase complex with the co-crystallized ligands used in the present study were retrieved from the Protein Data Bank (PDB code 3S41) and imported into HERMES 1.6 visualization interface implemented in Crystal Structure Theory and Applications the Cambridge Structural Database (CSD). Survey of the activator binding site was performed with the reference of amino acid residues of large GK domain as previously reported [9] [14]. Hydrogen atoms were then added to the protein for correct ionization and tautomeric states of amino acid residues. The molecular docking of 3S41 was carried out with GOLD 5.6 after conversion of its structure from 2D to 3D without minimization of protein energy.

Preparation of Ligands
The carboxamide ligand was the S41 co-crystallized activator of GK. One fragment of S41 structure was modified as labelled in Scheme 1 in order to design its 30 seleno, sulfo and oxoderivatives using Mercury 3.10 program implemented in CSD.    [19], using the 6 -31 G basis set.

Molecular Ligand-Protein Docking
The ligand-protein docking was carried out using the GOLD 5.6 program. Before docking GK with ligands, the reproduction of the experimental poses of S41 using 3S41 and 24 additional GK complexes with their co-crystallized ligands derived from PDB as shown in Table 2 was carried out to validate GOLD 5.6. A good reproduction was defined as one with a value of the root mean square deviation (RMSD) lower or equal to 2 Å [20] [21]. The default GOLD fitness function ChemPLP was used for the docking of GK-3S41 with derivatives by considering the flexibility of the ligands and specific amino acid residues in the activator site. The greater the ChemPLP fitness score, the better is the ligand binding affinity. Rescoring with Chemscore was performed to determine the energy of the ligand interaction with the GK site (ΔG bind ) according Scheme 2.
Free enthalpy of binding, where the ΔG 0 and S (bar) components are experimental and entropic terms, while ΔG HB , ΔG Met , ΔG Lip , ΔG Rot represent energy terms related to hydrogen-bonding, metal, lipophilic and rotation interactions, respectively [22].
Each simulation was performed ten times providing ten docked conformations until three of the ten poses were within 1.5 Å RMSD of each other. The ligand conformer pose with the lowest energy was considered as the binding conformation in the protein site. Figure 1 and Table 3, respectively, report the pose of the S41 in the activator site derived from PDB and the pose provided by GK docking with the GOLD 5.6 program.  [14]; (b) Pose of the co-crystallized S41 reproduced by GK docking with GOLD 5.6. Table 3. Experimental pose of the S41 ligand derived from PDB and its pose with geometric patterns of the specific interactions provided by docking using GOLD 5.6.  Table 3 show that the experimental pose of the co-crystallized ligand (S41) derived from PDB is well reproduced by docking with the GOLD 5.6 program. Indeed, the pattern of the ligand pose in Figure 1 matches well the activator site with its three hydrophobic pockets in form of the Y letter [6]. The   with RMSD values > 2 Å. Thus, most of RMSD values are equal or less than the limit value recommended for a good pose reproduction [20] [21] suggesting that GOLD 5.6 is a powerful docking program [18] [22].

Validation of GOLD 5.6 Potential for Ligand-Protein Docking
The experimental pose of co-crystallized activator from PDB (S41) and superposition of the poses derived from the molecular docking with GOLD 5.6 are illustrated in Figure 3.

Glucokinase Docking with Ligands
The areas of contacts between GK and ligands are illustrated in Figure 4.  Figure 4 reveals that hydrogen bonds and hydrophobic interactions were involved in the ligand binding as expected. This finding is consistent as the activator site is a combination of three hydrophobic pockets [6], yet it contains residues of polar amino acids, which are able to form hydrogen bonds with ligands [23]. Indeed, the first pocket contains Val62, Ile59, Val452 and Val255, which are hydrophobic, and the polar Arg63. The second pocket contains Pro66, Val65 and Tyr215, while the third pocket contains Met210, Met235 and Tyr214. Ligand areas allowing hydrogen-bonding and lipophilic interactions were evidenced by the molecular electrostatic potential surfaces (MEPs) depicted in Figure 5. It can be seen that MEPs of S41 and derivatives feature marked positive (in blue) and negative (in red) areas corresponding to hydrogen-bond donors or acceptors, but also large non-polar regions (in grey). Thus, lipophilic interactions are probably due to the contacts between hydrophobic amino acid residues present in the activator site of GK and non-polar areas of ligand molecules, while the hydrogen bonds can be ascribed to electrostatic interactions between hydrogen-bond donors and acceptors of polar amino acid (AA) residues and those of S41 and its derivatives. The affinity of the ligands to the binding site was then determined based on GOLD fitness scores, while the binding stability was obtained by calculating the energy of interaction [18]. Table 4 reports the values of GOLD fitness scores and energies of ligand-GK interaction. The values of GOLD fitness scores and free enthalpies of ligand binding suggest that the S41 activator and its derivatives were bonded with high affinity to GK. According to the fitness scores, 12 derivatives interact with GK more strongly than S41. All values of the binding enthalpies are negative suggesting that all ligand-GK complexes are thermodynamically stable. The binding enthalpies of the most stable complexes decrease in the sequence 3S> 5S > 1S >1O > 3O > 4O > 8S > 3Se. Values of energy component terms (ΔG HB and ΔG Lip ) show that the main forces stabilizing ligand-protein complexes are lipophilic interactions, enhanced by hydrogen bonding interactions. This is evidenced by the increase of ΔG Lip when an oxygen atom in the ligand is replaced by a chalcogen atom (S or Se). Such hydrogen-binding and lipophilic interactions should improve the activator selectivity and lead to least side effects, when they involve specific amino acid residues [6].
The predominance of the lipophilic interactions can be ascribed to hydrophobic amino acid residues present in the activator site. Indeed, it has been found that the Glu221, Met235, Ile211, Val455, Tyr215, Val91, Ala454, Leu451, Tyr214 and Val 62 amino acid residues were involved in the lipophilic interactions. Val62 and Tyr 215, respectively, interacted with ligands in the first and second pockets of the activator site, while Met235 and Tyr 214 were involved in hydrophobic effects in the third pocket as shown in Figure 4. This finding is in good agreement with the literature [6]. Table 5 reports distances and angles related to hydrogen bonds found in the binding site for the 6 best ranked ligands.
The geometric parameters were calculated according to Scheme 3. The value vdW (H) = 1.10 Å determined by Rowland et al. (1996) [24] was used for van der Waals radius of hydrogen atom, while the values vdW(O) = 1.52 Å for oxygen, vdW(N) = 1.55 Å for nitrogen and vdW(Se) = 1.90 Å for selenium were taken from Bondi et al. [25].
Scheme 3. Geometric parameters of the hydrogen-bonding: interaction distance (δ) and hydrogen-bond angle (θ). All values of the δ distance are less than the sum of van der Waals radii of the interacting atoms and the θ angles are greater than the 120˚ recommended limit [26]. This suggests that the hydrogen bonds were effective. Interestingly, all derivatives like the S41 ligand were able to form two specific hydrogen bonds with the Arg 63 amino acid residue, which is catalytically crucial for therapeutic properties of GK activators [6] [27].
It is noteworthy that beside hydrogen bonds and lipophilic interactions, ligands containing divalent S or Se atoms formed sigma-hole bonds, also known as chalcogen bonds as shown in Figure 6. The bond pattern in Figure 6 suggests that the chalcogen interaction were directed along the extensions of covalent bonds of the chalcogen atom as expected [28]. Table 6 reports the distances (δ) and angles (θ) of chalcogen bonds formed between the 3S, 3Se and 9Se ligands and lone pairs of atoms of the amino acid residues acting as bond acceptors in the binding site. As shown in Table 6, the mean value of 144˚ obtained for the bond angle (θ 1 ) was greater than 140˚ recommended limit for sigma-hole interactions [28]. The δ distances between interacting atoms are less than the sums of their respective Crystal Structure Theory and Applications

Conclusion
Investigations of the activation of glucokinase as a protein target in the treatment of diabetes have benefited from evident progress in drug design based on ligand-protein docking and quantitative structure-activity relationship (QSAR). In the present work, the docking of GK with seleno, sulfo and oxo derivatives of the S41 carboxamide activator revealed that all derivatives were bonded to the allosteric site of GK with relatively high affinity. According to ligand binding energies, the main forces stabilizing the complexes were lipophilic interactions enhanced by hydrogen bonds, but also chalcogen bonds in the case of seleno and sulfo derivatives. The increase of thermodynamic stability of these complexes compared with 3S41 can be ascribed to enhancement of the lipophilic character of S41 derivatives. Interestingly, twelve of the S41 derivatives showed stronger binding affinity than the co-crystallized ligand, while maintaining the two hydrogen bonds with crucial Arg63 amino acid residue. Therefore, this study provides an additional support to recent advances in use of the molecular docking as a powerful tool for design and discovery of activators with improved pharmacological properties. The study of the twelve best ranked derivatives should be extended to pharmacophore investigations on their therapeutic properties. The prediction of their activator activities through molecular dynamic simulation and QSAR studies based on ligand chemical and quantum descriptors should also be useful for experimental and clinical trials.