Computational Study on New Natural Compound Inhibitors of Pyruvate Dehydrogenase Kinases

Pyruvate dehydrogenase kinases (PDKs) are key enzymes in glucose metabolism, negatively regulating pyruvate dehyrogenase complex (PDC) activity through phosphorylation. Inhibiting PDKs could upregulate PDC activity and drive cells into more aerobic metabolism. Therefore, PDKs are potential targets for metabolism related diseases, such as cancers and diabetes. In this study, a series of computer-aided virtual screening techniques were utilized to discover potential inhibitors of PDKs. Structure-based screening using Libdock was carried out following by ADME (adsorption, distribution, metabolism, excretion) and toxicity prediction. Molecular docking was used to analyze the binding mechanism between these compounds and PDKs. Molecular dynamic simulation was utilized to confirm the stability of potential compound binding. From the computational results, two novel natural coumarins compounds (ZINC12296427 and ZINC12389251) from the ZINC database were found binding to PDKs with favorable interaction energy and predicted to be non-toxic. Our study provide valuable information of PDK-coumarins binding mechanisms in PDK inhibitor-based drug discovery.


Introduction
Pyruvate dehydrogenase kinases (PDKs) are key enzymes in glucose metabolism since they negatively regulate the activity of pyruvate dehydrogenase complex (PDC) by phosphorylation [1,2]. PDC is an important gatekeeper enzyme that catalyzes the oxidative decarboxylation of pyruvate to produce acetyl CoA [3]. The mammalian PDC is a 9.5 million Dalton multi-enzyme complex consisting of pyruvate dehydrogenase (E1), dihydrolipoyltransacetylase (E2), dihydrolipoamide dehydrogenase (E3) and the E3-bindingprotein (E3BP). PDKs are non-covalently bound to the lipoyl domain (L2) of the E2 subunit of PDC. They phosphorylate the pyruvate dehydrogenase (PDH)-E1α and inactivate PDC [4][5][6][7]. In mammalians, there are four isoenzymes of PDKs (PDK 1 to 4) with different binding affinity, phosphorylation site specificity and tissue distribution [8,9]. The binding affinities of four isoenzymes to L2 are in the following order: PDK3 > PDK1 « PDK2 > PDK4. As for the phosphorylation sites (S264, S271 and S203), only PDK1 is reported to phosphorylate all three sites, while other isoenzymes can only phosphorylate S264 and S271 with different rates. Four isoenzymes have different tissue distribution: PDK1 is present preferentially in the pancreatic islets, heart and skeletal muscles; PDK2 is expressed universally in all tissues; PDK3 is abundant in the kidney, brain and testes; PDK4 is more easily detected in the pancreatic islets, heart, kidney and skeletal muscle.

Virtual Screening of Natural Products Database against Pyruvate Dehydrogenase Kinase 1 (PDK1)
Lipoamide-binding pocket is an important regulatory site of PDKs, since small molecules binding to this pocket can prevent PDKs from binding to PDC but with minimal effect on the function of the enzyme complex. To identify new compounds that could potentially inhibit PDK1 through binding to the lipoamide-binding pocket, a virtual screening was carried out using Libdock [35] module of Discovery Studio 3.5 (DS3.5, Accelrys, Inc., San Diego, CA, USA). 99,932 purchasable natural product molecules were taken from the ZINC database, which is a free database of commercially-available compounds provided by the Irwin and Shoichet Laboratories in the Department of Pharmaceutical Chemistry at the University of California, San Francisco (UCSF). The target used in this study was the 3D structure of PDK1 (PDB ID: 2Q8G) [25]. The well known PDK inhibitor AZD7545, which can inhibit PDKs activity except PDK4 in vitro and in vivo [24,36], was chosen as reference compound. AZD7545 inhibited PDK1 through the trifluoromethylpropanamide end that inserted into the lipoamide-binding pocket of PDK1, as revealed by the crystal structure of human PDK1-AZD7545 complex. The blocking of the lipoamide-binding pocket resulted in inhibition of PDKs activities by aborting kinase binding to the PDC scaffold [25,37,38]. After the screening, 2354 compounds were found have higher Libdock scores than AZD7545 (Libdock score: 117.276). The top 20 ranked compounds are listed in Table 1.
Safety is an important aspect of drug research. To examine the safety of the compounds, different toxicity such as Ames mutagenicity (AMES), rodent carcinogenicity (based on the U.S. National Toxicology Program (NTP) dataset) and developmental toxicity potential (DTP) properties of compounds and AZD7545 were predicted using TOPKAT module of DS (Table 3). The results showed that, all compounds were predicted to be non-mutagen. Eleven compounds were predicted to be non-carcinogen and seven compounds with no developmental toxicity potential. The reference AZD7545 was predicted with developmental toxicity potential. Synthesizing the above results, ZINC12296427 and ZINC12389251 are not CYP2D6 inhibitors, with no hepatotoxicity. Moreover, they are predicted with no Ames mutagenicity, rodent carcinogenicity and developmental toxicity potential. Therefore, ZINC12296427 and ZINC12389251 were predicted safe drug candidates and selected for further research (Figure 1).  Safety is an important aspect of drug research. To examine the safety of the compounds, different toxicity such as Ames mutagenicity (AMES), rodent carcinogenicity (based on the U.S. National Toxicology Program (NTP) dataset) and developmental toxicity potential (DTP) properties of compounds and AZD7545 were predicted using TOPKAT module of DS (Table 3). The results showed that, all compounds were predicted to be non-mutagen. Eleven compounds were predicted to be non-carcinogen and seven compounds with no developmental toxicity potential. The reference AZD7545 was predicted with developmental toxicity potential. Synthesizing the above results, ZINC12296427 and ZINC12389251 are not CYP2D6 inhibitors, with no hepatotoxicity. Moreover, they are predicted with no Ames mutagenicity, rodent carcinogenicity and developmental toxicity potential. Therefore, ZINC12296427 and ZINC12389251 were predicted safe drug candidates and selected for further research (Figure 1).

Ligand Binding Analysis
In order to study ligand binding mechanisms, ZINC12296427 and ZINC12389251 were docked into the 3D structure of PDK1 by using the CDOCKER module of DS [39]. The CDOCKER interaction energy as an estimation of molecular complex binding affinity was used in this study. As shown in Table 4, the CDOCKER interaction energy of ZINC12296427 is lower than the reference AZD7545, indicating that ZINC12296427 may have a higher binding affinity with PDK1. The energy of the other compound ZINC12389251 is close to AZD7545. We analyzed the hydrogen bonds and Pi-Pi interactions formed by PDK1 and these compounds (Figures 2 and 3, Tables 5 and 6). The results showed that, both ZINC12296427 and ZINC12389251 formed one pair of hydrogen bond with PDK1, by the carbonyl O of compounds and the amide -NH of GLN197 of PDK1. No hydrogen bonds were formed between AZD7545 and PDK1. Each of the three compounds formed one pair Pi-Pi interaction with PDK1, by the centroids of the benzene ring of PHE65 and the centroids of the benzene ring of coumarin groups of ZINC12296427 and ZINC12389251 or chlorinated benzene ring of AZD7545, respectively. The additional hydrogen bond formed between ZINC12296427, ZINC12389251 and the lipoamide-binding pocket of PDK1might contribute to the lower interaction energy of the complexes. These results indicated that, ZINC12296427 and ZINC12389251 probably bind to PDK1 with similar or even better affinity than AZD7545.

Ligand Binding Analysis
In order to study ligand binding mechanisms, ZINC12296427 and ZINC12389251 were docked into the 3D structure of PDK1 by using the CDOCKER module of DS [39]. The CDOCKER interaction energy as an estimation of molecular complex binding affinity was used in this study. As shown in Table 4, the CDOCKER interaction energy of ZINC12296427 is lower than the reference AZD7545, indicating that ZINC12296427 may have a higher binding affinity with PDK1. The energy of the other compound ZINC12389251 is close to AZD7545. We analyzed the hydrogen bonds and Pi-Pi interactions formed by PDK1 and these compounds (Figures 2 and 3 Tables 5 and 6). The results showed that, both ZINC12296427 and ZINC12389251 formed one pair of hydrogen bond with PDK1, by the carbonyl O of compounds and the amide -NH of GLN197 of PDK1. No hydrogen bonds were formed between AZD7545 and PDK1. Each of the three compounds formed one pair Pi-Pi interaction with PDK1, by the centroids of the benzene ring of PHE65 and the centroids of the benzene ring of coumarin groups of ZINC12296427 and ZINC12389251 or chlorinated benzene ring of AZD7545, respectively. The additional hydrogen bond formed between ZINC12296427, ZINC12389251 and the lipoamide-binding pocket of PDK1might contribute to the lower interaction energy of the complexes. These results indicated that, ZINC12296427 and ZINC12389251 probably bind to PDK1 with similar or even better affinity than AZD7545.

Ligand Binding Analysis
In order to study ligand binding mechanisms, ZINC12296427 and ZINC12389251 were docked into the 3D structure of PDK1 by using the CDOCKER module of DS [39]. The CDOCKER interaction energy as an estimation of molecular complex binding affinity was used in this study. As shown in Table 4, the CDOCKER interaction energy of ZINC12296427 is lower than the reference AZD7545, indicating that ZINC12296427 may have a higher binding affinity with PDK1. The energy of the other compound ZINC12389251 is close to AZD7545. We analyzed the hydrogen bonds and Pi-Pi interactions formed by PDK1 and these compounds (Figures 2 and 3, Tables 5 and 6). The results showed that, both ZINC12296427 and ZINC12389251 formed one pair of hydrogen bond with PDK1, by the carbonyl O of compounds and the amide -NH of GLN197 of PDK1. No hydrogen bonds were formed between AZD7545 and PDK1. Each of the three compounds formed one pair Pi-Pi interaction with PDK1, by the centroids of the benzene ring of PHE65 and the centroids of the benzene ring of coumarin groups of ZINC12296427 and ZINC12389251 or chlorinated benzene ring of AZD7545, respectively. The additional hydrogen bond formed between ZINC12296427, ZINC12389251 and the lipoamide-binding pocket of PDK1might contribute to the lower interaction energy of the complexes. These results indicated that, ZINC12296427 and ZINC12389251 probably bind to PDK1 with similar or even better affinity than AZD7545.      To confirm the docking analysis performed with CDOCKER, the docking results were crosschecked using AutoDock [40]. All docking conformations were visualized using DS3.5 so as to ensure the ligands were docked into the defined binding pocket. The superimposed docked structures are shown in Figure 4. For ZINC12296427 and AZD7545, the docked poses generated by AutoDock were similar to the poses generated by CDOCKER, with RMSD (Root Mean Square Deviation) values of 1.5 and 0.3 Å, respectively. ZINC12296427 had a lower binding energy (´8.8 kcal/mol) than AZD7545 (´6.3 kcal/mol) ( Table 7), in accordance with the results of CDOCKER. There were some differences between the docked poses of ZINC12389251 generated by AutoDock and CDOCKER. In the PDK1-ZINC12389251 complex predicted by AutoDock, there was one more pair of Pi-Pi interaction compared to the results of CDOCKER, formed by the centroids of the benzene ring of PHE65 and the coumarin group of ZINC12389251.This extra pair of Pi-Pi interaction might contribute to the affinity of ZINC12389251 with PDK1, making ZINC12389251 have the lowest binding energy (´9.3 kcal/mol) in three compounds (Table 7). Overall, both ZINC12296427 and ZINC12389251 bind to lipoamide binding pocket of PDK1 with lower energy than AZD7545 predicted by the two methods.
There were some differences between the docked poses of ZINC12389251 generated by AutoDock and CDOCKER. In the PDK1-ZINC12389251 complex predicted by AutoDock, there was one more pair of Pi-Pi interaction compared to the results of CDOCKER, formed by the centroids of the benzene ring of PHE65 and the coumarin group of ZINC12389251.This extra pair of Pi-Pi interaction might contribute to the affinity of ZINC12389251 with PDK1, making ZINC12389251 have the lowest binding energy (−9.3 kcal/mol) in three compounds (Table 7). Overall, both ZINC12296427 and ZINC12389251 bind to lipoamide binding pocket of PDK1 with lower energy than AZD7545 predicted by the two methods.

Molecular Dynamics Simulation
To evaluate the stabilities of PDK1-ligand complexes under dynamic conditions, we conducted the molecular dynamics (MD) simulation using DS. The initial conformations were acquired from the molecular docking experiments by CDOCKER. The RMSD curves of the receptor structures from each complex and the potential energy profiles of each complex are shown in Figure 5. The trajectories of complexes reached equilibrium after 10 ns, RMSD and potential energy of the complexes gets

Molecular Dynamics Simulation
To evaluate the stabilities of PDK1-ligand complexes under dynamic conditions, we conducted the molecular dynamics (MD) simulation using DS. The initial conformations were acquired from the molecular docking experiments by CDOCKER. The RMSD curves of the receptor structures from each complex and the potential energy profiles of each complex are shown in Figure 5. The trajectories of complexes reached equilibrium after 10 ns, RMSD and potential energy of the complexes gets stabilized with the time. The H-bond distances formed between ZINC12296427, ZINC12389251 and PDK1 were within a range around 3.0 and 2.3 Å. In addition, both ZINC12296427 and ZINC12389251 formed two pairs of hydrogen bonds with the water molecules after molecular dynamics simulation (TIP21914:H2-ZINC12296427:O21; ZINC12296427:H28-TIP13772:OH2; TIP12999:H2-ZINC12389251:O18; ZINC12389251:H45-TIP17204:OH2). These hydrogen bonds might contribute to the stability of the complexes. Taking all the evaluation indexes into consideration, these two compounds might interact with PDK1 steadily and have potential negative modulatory effects on PDK1.

Molecular Docking with Other PDKs
ZINC12296427 and ZINC12389251 were docked into the 3D structures of PDK2, PDK3 and PDK4 to evaluate if they could inhibit other PDKs. CDOCKER and AutoDock results showed that, except for PDK4, both compounds could bind to the lipoamide binding pockets of PDK2 and PDK3 with strong interaction energy (Tables 4 and 7), and formed hydrogen bond (Table 5) and Pi-Pi interactions ( Table 6) with receptors. These results indicated that ZINC12296427 and ZINC12389251 might be broad spectrum inhibitors of PDKs.

Discussion
To phosphorylate the pyruvate dehydrogenase (E1 of PDC), PDKs need to interact with the L2 domain of dihydrolipoyltransacetylase (E2 of PDC) through their lipoamide binding pockets. Inhibitors binding to the lipoamide binding pockets can uncouple the link between PDKs and PDC, and then upregulate PDC activity indirectly. In this study, we identified two novel natural compounds, ZINC12296427 and ZINC12389251, from ZINC database to effectively inhibit PDKs through virtual screening technique. Molecular docking study showed that the two novel compounds could bind to the lipoamide-binding pocket of PDKs through forming of hydrogen bonds and Pi-Pi interactions with receptors. Moreover, these two compounds were predicted with no hepatotoxicity, Ames mutagenicity, rodent carcinogenicity and developmental toxicity. Thus, our findings provided potential mechanism of two novel inhibitors of PDKs for clinical drug design. Since the compounds were virtually screened and their inhibitory activities have not been reported, experiments such as IC50 and EC50 measurements will be carried out in our further study to examine their bioactivities.
Our results of the binding mechanism study also showed that the functional groups of these two compounds are their coumarins group. Coincidentally, recent studies have reported the antitumor activities of coumarins and their derivatives [41]. Our results suggest that upregulation of mitochondria aerobic metabolism of cancer cells through inhibition of PDKs might contribute to the anticancer activities of this kind of compound. This is also probably one possible mechanism for the anticancer activity of coumarones, which have similar structure to coumarins. Further studies are needed to determine the exact mechanism. Coumarins have broad pharmacological activities, such

Molecular Docking with Other PDKs
ZINC12296427 and ZINC12389251 were docked into the 3D structures of PDK2, PDK3 and PDK4 to evaluate if they could inhibit other PDKs. CDOCKER and AutoDock results showed that, except for PDK4, both compounds could bind to the lipoamide binding pockets of PDK2 and PDK3 with strong interaction energy (Tables 4 and 7), and formed hydrogen bond (Table 5) and Pi-Pi interactions ( Table 6) with receptors. These results indicated that ZINC12296427 and ZINC12389251 might be broad spectrum inhibitors of PDKs.

Discussion
To phosphorylate the pyruvate dehydrogenase (E1 of PDC), PDKs need to interact with the L2 domain of dihydrolipoyltransacetylase (E2 of PDC) through their lipoamide binding pockets. Inhibitors binding to the lipoamide binding pockets can uncouple the link between PDKs and PDC, and then upregulate PDC activity indirectly. In this study, we identified two novel natural compounds, ZINC12296427 and ZINC12389251, from ZINC database to effectively inhibit PDKs through virtual screening technique. Molecular docking study showed that the two novel compounds could bind to the lipoamide-binding pocket of PDKs through forming of hydrogen bonds and Pi-Pi interactions with receptors. Moreover, these two compounds were predicted with no hepatotoxicity, Ames mutagenicity, rodent carcinogenicity and developmental toxicity. Thus, our findings provided potential mechanism of two novel inhibitors of PDKs for clinical drug design. Since the compounds were virtually screened and their inhibitory activities have not been reported, experiments such as IC 50 and EC 50 measurements will be carried out in our further study to examine their bioactivities.
Our results of the binding mechanism study also showed that the functional groups of these two compounds are their coumarins group. Coincidentally, recent studies have reported the antitumor activities of coumarins and their derivatives [41]. Our results suggest that upregulation of mitochondria aerobic metabolism of cancer cells through inhibition of PDKs might contribute to the anticancer activities of this kind of compound. This is also probably one possible mechanism for the anticancer activity of coumarones, which have similar structure to coumarins. Further studies are needed to determine the exact mechanism. Coumarins have broad pharmacological activities, such as anti-coagulant, anti-viral, anti-inflammatory and anti-microbial effects [42]. Therefore, development of novel bioactivities of natural and synthetic coumarins based on their PDKs inhibition activities is of great significance for the treatment of metabolic diseases.

Docking Software and Ligand Library
Libdock, and ADMET modules of Discovery Studio 3.5 software (DS3.5, Accelrys, Inc.) were used for virtual screening and ADMET analyzing of inhibitors. CDOCKER and AutoDock were used for docking study. A Natural Products database (NP) in the ZINC database was employed to screen PDKs inhibitors.

Structure-Based Virtual Screening Using Libdock
Libdock is a rigid-based docking program. It calculates hotspots for the protein using a grid placed into the binding site and using polar and apolar probes (San Diego, CA, USA). Then the hot spots are further used to align the ligands to form favourable interaction. The Smart Minimiser algorithm and CHARMm force field (Cambridge, MA, USA) were used for ligands minimization. After minimized, all the ligand poses are ranked based on the ligands score. The 1.9 Å crystal structure of human pyruvate dehydrogenase kinase 1 (PDK1) in complex with AZD7545 (PDB ID: 2Q8G) were downloaded from protein data bank (PDB) and imported to the working environment of Libdock. The protein was prepared by removing crystal water and other hetero atoms (except AZD7545), followed by addition of hydrogen, protonation, ionization and energy minimization. The CHARMm [43] force field and the Smart Minimiser algorithm were applied for energy minimization. The minimization performed 2000 steps with an RMS gradient tolerance of 0.1, and the final RMS gradient was 0.09463. The prepared protein was used to define the binding site from the "Edit binding site" option on the receptor-ligand interaction tool bar. Using the bound ligands (AZD7545) binding positions, the active sites for docking were generated. Virtual screening was carried out by docking all the prepared ligands at the defined active site using Libdock. Based on the Libdock score, all the docked poses were ranked and grouped by name. All compounds were ranked according to their Libdock score.

Molecular Docking
CDOCKER module of Discovery Studio was used for molecular docking study. CDOCKER is an implementation of a CHARMm based docking tool. The receptor is held rigid while the ligands are allowed to flex during the docking process. For each complex pose, the CHARMm energy (interaction energy plus ligand strain) and the interaction energy, which indicate ligand binding affinity, are calculated. Crystal structure of PDK1 (PDB ID: 2Q8G, 1.9 Å), PDK2 (PDB ID: 2BU6, 2.4 Å) [44], PDK3 (PDB ID: 2Q8I, 2.6 Å) [25] and PDK4 (PDB ID: 2ZDX, 2.5 Å) [26] were obtained from the protein data bank. The crystal water molecules are generally removed in rigid and semi-flexible docking process [45,46], since the fixed water molecules might affect the formation of receptor-ligand complex. The water molecules were removed and hydrogen atoms were added to the protein. The 3D structures of ZINC12296427 and ZINC12389251 were obtained from ZINC database (no experimental crystal structures available in crystallographic databases). The CHARMm forcefield was used for receptors and ligands. The binding site spheres of PDK1 and PDK2 were defined as the regions that come within radius 15 Å from the geometric centroid of the ligands AZD7545, AZ12, respectively. The binding site spheres of PDK3 and PDK4 were defined by structural alignment with PDK1. During the docking process, the ligands were allowed to bind to the residues within the binding site spheres. After being extracted from the binding site, the initial compound AZD7545 was re-docked into the crystal structure of PDK1. The RMSD between the docked pose and the crystal structure of the complex was 0.6 Å, indicating the CDOCKER module was highly reliable for reproducing the experimentally observed binding mode of PDKs. The structures of identified hits were prepared and docked into the lipoamide-binding pocket of PDKs. Different poses for each test molecule were generated and analyzed on the basis of CDOCKER interaction energy.
AutoDock 4.2 (The Scripps Research Institute, (La Jolla, CA, USA) was used to crosscheck the docking results of CDOCKER. The AutoDockTools 1.5.6 was used to prepare the protein and ligands for docking procedure. All solvent molecules, water molecules, and the cocrystallized ligand were removed from the structure; Kollman charges and polar hydrogens were added. AutoGrid (La Jolla, CA, USA) was used to generate the grid maps. Each grid was centered at the lipoamid binding site of the PDKs. The grid dimensions were 45 points in each dimension separated by 0.375 Å. The files were generated as PDBQT format. For all ligands, random starting positions, orientations and torsions were used. PDBQT file of the ligands was generated with all the default values accepted. The Genetic Algorithm was used with 2,500,000 energy evaluations and a population of 300 individuals; 100 runs were carried out. Following docking, the results were clustered into groups with RMSD <2.0 Å. The ranking of the clusters was performed based on the lowest energy representative of each cluster. Visualization and analysis of interactions between protein and ligand were performed using DS3.5.

Molecular Dynamics Simulation
The best binding conformations of the PDK1-inhibitor complexes among the poses predicted by the molecular docking program were selected and used as MD simulation starting points. The ligand-receptor complex was put into an orthorhombic box and solvated with an Explicit Periodic Boundary solvation water model. The size of the box was 88ˆ75ˆ72 Å, containing 13,802 water molecules. In order to simulate the physiological environment, sodium chloride were added to the system with the ionic strength of 0.145. Then, the system was subjected to the CHARMm force-field and relaxed by energy minimization (1000 steps of steepest descent and 1000 steps of conjugated gradient), with the final RMS gradient of 0.08326. The system was slowly driven from an initial temperature of 50 K to the target temperature of 300 K for 2 ns and equilibration simulations were run for 1 ns. The MD simulations (production) were performed for 40 ns with time step of 2 fs. The simulation was performed with the NPT (normal pressure and temperature) system at a constant temperature of 300 K and the results were saved at a frequency of 0.02 ns. The Particle Mesh Ewald (PME) algorithm was used to calculate long range electrostatics, and the linear constraint solver (LINCS) algorithm was adapted to fix all bonds involving hydrogen. The initial complex was set as a reference. The MD trajectory was determined for structural properties, root mean-square deviation (RMSD), and potential energy by using the Discovery Studio 3.5 analyze trajectory protocol (San Diego, CA, USA).

ADME and Toxicity Prediction
The ADMET and TOPKAT modules of DS were employed to calculate the ADMET properties of the compounds, such as their aqueous solubility, blood-brain barrier (BBB) penetration, cytochrome P450 2D6 (CYP2D6) inhibition, hepatotoxicity, human intestinal absorption, plasma protein binding (PPB) level, rodent carcinogenicity, Ames mutagenicity and developmental toxicity potential.

Conclusions
The present study elucidates the possible inhibitory activity of two natural compounds ZINC12296427 and ZINC12389251 target PDKs. The ADMET properties of these compounds were predicted to be non-carcinogenic and non-toxic, indicating safe drug candidates. The results of this study not only demonstrate the probable binding mode of these compounds with PDKs, but also encourage further investigations of the effect of coumarins on metabolism through inhibition of PDKs.