In Silico Design for Adenosine Monophosphate-Activated Protein Kinase Agonist from Traditional Chinese Medicine for Treatment of Metabolic Syndromes

Adenosine monophosphate-activated protein kinase (AMPK) acts as a master mediator of metabolic homeostasis. It is considered as a significant millstone to treat metabolic syndromes including obesity, diabetes, and fatty liver. It can sense cellular energy or nutrient status by switching on the catabolic pathways. Investigation of AMPK has new findings recently. AMPK can inhibit cell growth by the way of autophagy. Thus AMPK has become a hot target for small molecular drug design of tumor inhibition. Activation of AMPK must undergo certain extent change of the structure. Through the methods of structure-based virtual screening and molecular dynamics simulation, we attempted to find out appropriate small compounds from the world's largest TCM Database@Taiwan that had the ability to activate the function of AMPK. Finally, we found that two TCM compounds, eugenyl_beta-D-glucopyranoside and 6-O-cinnamoyl-D-glucopyranose, had the qualification to be AMPK agonist.


Introduction
One study has found that the cell "starvation" signal transduction pathway reveals the process of cell "starvation" signal transduction mechanisms. This finding is considered as a significant millstone to treat metabolic syndromes including obesity, diabetes, and fatty liver [1,2]. Adenosine monophosphate-activated protein kinase (AMPK) was always a hot issue in the recent 10 years because it plays an essential role in sensing cellular energy or nutrient status [3,4]. AMPK exists in several organs including skeletal muscles, brain, and liver [5]. For example, when skeletal muscles have exercised for a period of time, activated AMPK helps the cells adapt to energy challenge by increasing glucose uptake [6]. AMPK acts as a master mediator of metabolic homeostasis [7]. The effect of AMPK activation is switching on the catabolic pathways consisting of cellular glucose uptake, glycolysis, fatty acid oxidation, inhibition of triglyceride, cholesterol, and protein synthesis [8,9]. The upstream signal molecule, such as adiponectin, stimulates AMPK to utilize glucose and fatty acid [10]. The activated AMPK in the hypothalamus increases the appetite or desire for food intake [11,12]. Adiponectin stimulates AMPK to protect the heart against myocardial ischemia [13,14]. Another upstream signal molecule, leptin, activates AMPK to stimulate fatty acid oxidation [15]. AMPK in the hypothalamus mediates thyroid hormone to regulate energy balance, too [16].
AMPK consists of three subunits, , , and [17]. Each of the subunits has a specific structure and function [18,19]. The following mechanism explains how the three subunits work. When the cells need energy, adenosine triphosphate (ATP) converts to adenosine diphosphate (ADP) and AMP by releasing energy simultaneously. At the normal condition, excess AMP binds to the subunit of AMPK [20]. After binding of subunit exposes the active site, Thr172, on the catalytic subunit of AMPK [21]. For all of the three subunits, subunit is the most important. Activation of AMPK follows conformational change of subunit and phosphorylation at Thr172 [22]. The phosphorylation is induced by the upstream signal molecule such as liver kinase B1 (LKB1) [   subunit is located between and subunits and is associated with the function of glycogen sensor [24]. ADP can also bind to subunit and protect AMPK from dephosphorylation but cannot cause conformational change [25]. Physiological study has approved that blood glucose is partly regulated by AMPK [26]. Activating AMPK by its agonist decreases insulin resistance induced by obesity [27]. Drug design for AMPK activator provides the hope for treating type 2 diabetes mellitus [28]. AMPK has become the drug target for managing diabetes and metabolic syndrome [29]. Metformin is the classic antidiabetes medicine involved in activation of AMPK [30,31]. Investigation of AMPK has new findings recently. Activated AMPK has the function of anti-inflammation by a mimic state of pseudostarvation [32]. AMPK can regulate cell growth via signal integration [33]. The LKB1-AMPK pathway can mediate cellular autophagy or apoptosis by phosphorylating p27kip1 [34]. Thus the LKB1-AMPK pathway has been approved in the field of tumor suppression [35]. AMPK inducing phosphorylation of Unc-51-like kinase 1 (ULK1) can cause autophagy in an energetic exhaustion status. ULK1 is a serine/threonine-protein kinase [36,37].
Due to progress of modern technology, the binding phenomena of protein dynamics motion and structure changing can be analyzed by computational simulation [38,39]. In silico investigation of biology or computational systems biology helps us to explore the protein-protein or protein-molecule interaction [40,41]. These methods make it possible to participate in computer-aided drug design (CADD) [42,43]. CADD techniques save our time to select appropriate drug compound rapidly compared with traditional one-by-one biochemistry [44,45]. CADD procedures consist of virtual screening, validations, and analysis [46,47]. Virtual screening and analysis employ docking and molecular dynamics (MD) simulation [48][49][50]. MD can predict how long the compound needs to form stable complex structure with target protein [51]. A series of statistics or score systems facilitate docking and MD accuracy [52]. Best candidates from virtual screening and MD can be selected as potential therapeutic drugs [53].
With progress in medical technology, many diseases can be resolved nowadays [54][55][56]. When we know how the diseases happen, we can handle them in smart ways [57][58][59]. Because of coordinating cell growth and autophagy, AMPK has become an emerging target for drug design of tumor inhibition [60]. Design of small molecular AMPK agonist is possible [61]. Small molecular drug design by CADD has been applied in systems biology extensively [62][63][64]. Traditional Chinese medicine (TCM) has therapeutic effect on many diseases [65][66][67][68]. In this study, we attempted to find out appropriate small compounds from the world's largest TCM Database@Taiwan that had the ability to activate the function of AMPK [69].

Compound Database.
To identify potential AMPK activators from TCM, we downloaded all small molecules from

Structure-Based Virtual
Screening. The ligands from TCM Database@Taiwan and the control ligand, adenosine monophosphate (AMP), were prepared for specified techniques. Ligand docking was performed by the LigandFit module in DS 2.5. The force field of Chemistry at HARvard Molecular Mechanics (CHARMm) was utilized to minimize all docking poses [70]. Absorption, distribution, metabolism, excretion, and toxicity (ADMET) were applied for all TCM compounds [71]. The scoring systems such as Dock score, piecewise linear potentials (PLP), and potential of mean force (PMF) were calculated by the LigandFit module in DS 2.5 [72,73].

Disorder Prediction.
We drew disorder disposition to exclude disordered residues by the program of PONDR-FIT in the DisProt website [74,75].

Molecular Dynamics (MD)
Simulation. The package of GROMACS (GROningen MAchine for Chemical Simulations) was used for MD simulation. We employed SwissParam to determine topology and parameters of small compounds for GROMACS simulation. The cytoplasmic condition was set with transferable intermolecular potential 3P (TIP3P) water at 0.9% sodium chloride concentration. After docking, selected protein-ligand complexes were conducted under the following phases: minimization, heating, equilibration, and production. The minimization protocol included 500 steps of steepest descent and 500 steps of conjugated gradient. The heating time was 50 ps from 50 K to 310 K. The equilibration time was 150 ps at 310 K. The production time was 5000 ps with constant temperature dynamics method. The time for temperature decay was 0.4 ps. We utilized the trajectory analysis to illustrate root mean square deviation (RMSD), Gyrate, solvent accessible surface (SAS), root mean square fluctuation (RMSF), total energy, database of secondary structure assignment (DSSP), matrices of smallest distance of residues for individual ligands, and protein during MD. We calculated the formation and distance of hydrogen bond (Hbond), too. Best distance of H-bond was set at 0.3 nm.

Ligand Pathway.
To analyze the ligand pathway, we employed the software of LigandPath module to illustrate the possible pathway of each ligand. A surface probe and minimum clearance were set at 6Å and 3Å, respectively.

Disorder Prediction.
Main residues (Val18 to Asp151) of AMPK-modeled structure for the 2 candidates and the control were not in the disordered area, so there was no influence to the shape of the binding site ( Figure 6).

Molecular Dynamics (MD)
Simulation. MD trajectories generated by GROMACS were illustrated. We utilized root mean square deviation (RMSD) to show the deviation degree of each ligand or protein from the beginning to the end of MD. When the ligand formed a complex with the protein (AMPK), eugenyl beta-D-glucopyranoside, 6-O-cinnamoyl-D-glucopyranose, and the control (AMP) had first larger deviations at 1600 ps, 1450 ps, and 2000 ps, respectively. AMP had a larger average deviation than the candidates (Figure 7(a)). We utilized Gyrate to measure the average distance of the atoms to the center of each mass. In other words, Gyrate showed the compact degree of each ligand or protein. 6-O-Cinnamoyl-D-glucopyranose had the largest ligand Gyrate values, but eugenyl beta-D-glucopyranoside had the largest protein Gyrate values. All the proteins had a trend to compact status at the end of MD (Figure 7(b)). We utilized solvent accessible surface (SAS) to measure  MD. All the candidates and the control had similar findings. Among the components of important second degree structure, the ratio of -helix was higher than -sheet, and -sheet was higher than bend. The ratio of -helix increased slightly, but the ratio of bend decreased comparatively ( Figure 10). We utilized matrices of the smallest distance to find any variation of residue distance when the protein complexed with the ligand. There was not any significant difference compared to the candidates with the control (Figure 11). Table 2 showed occupancy of H-bond between eugenyl beta-D-glucopyranoside, 6-O-cinnamoyl-Dglucopyranose, and AMP with AMPK protein during MD. The formation of H-bond between eugenyl beta-Dglucopyranoside with Asn138 and Asp151, 6-O-cinnamoyl-D-glucopyranose with Val18, Asp133, Lys135, and Asp151, and AMP with Gly19, Val90, and Glu94 was consistent with the results of docking. AMP also formed H-bond with Asp151, so we could conclude that Asp151 was key residue for top 2 TCM candidates and the control bound with AMPK.
We showed the distance of H-bonds between eugenyl beta-D-glucopyranoside and essential amino acids of AMPK. The O13 of eugenyl beta-D-glucopyranoside

Ligand
Pathway. 3D simulation of ligand pathway bound with AMPK helped us understand the process of combination. Eugenyl beta-D-glucopyranoside, 6-O-cinnamoyl-D-glucopyranose, and AMP had different pathways bound with AMPK protein (Figure 15).

Homology Modeling.
Because the catalytic subunit is the most important for all the three subunits of AMPK, we selected human AMPK (Q13131, subunit) and rat templates (2Y94, subunit) for homology modeling. The high percentage of identity and similarity of sequence alignment between Q13131 and 2Y94 indicated that the sequence alignment was reliable. The high percentage of residues in the favored and allowed area indicated that the AMPK-modeled structure was reasonable. The Verify scores of most residues were positive which indicated that the AMPK-modeled structure was reliable. We estimated that the negative values from residue 290 to residue 360 might be attributed to loss of reasonable structure from residue 311 to residue 341 illustrated in Figure 1.

Structure-Based Virtual Screening.
Whether absorption level of ADMET, Dock scores, PLP1, PLP2 or PMF, all the top 7 TCM compounds were better than the control (AMP). Absorption level of ADMET from 0 to 3 means good to very low absorption. We selected the first 2 compounds, eugenyl beta-D-glucopyranoside and 6-O-cinnamoyl-Dglucopyranose, as candidates due to their good absorption level. Both eugenyl beta-D-glucopyranoside and 6-Ocinnamoyl-D-glucopyranose formed H-bond with Asp151. Both eugenyl beta-D-glucopyranoside and AMP formed H-bond with Glu94. The main residues bound for the candidates and the control were located from Val18 to Asp151. The results mean that loss of reasonable structure from residue 311 to residue 341 did not disturb us for adopting the AMPK-modeled structure.  Figure 14: Distance of hydrogen bonds between AMP and essential amino acids of AMPK. All the compounds had larger fluctuations at the range between residue 320 and residue 400 in RMSF. The larger fluctuations were not reliable due to loss of reasonable structure from residue 311 to residue 341. Based on the similar graph in RMSF and total energy, we could speculate that eugenyl beta-D-glucopyranoside, 6-O-cinnamoyl-D-glucopyranose, and the control (AMP) had the same ability to activate AMPK. Our hypothesis was further approved by DSSP and smallest distance matrices. The structure of AMPK had similar change when it complexed with the 2 candidates and the control. Eugenyl beta-D-glucopyranoside, 6-O-cinnamoyl-D-glucopyranose, and the control had the ability to induce AMPK conformational change.
According to occupancy of H-bond between eugenyl beta-D-glucopyranoside, 6-O-cinnamoyl-Dglucopyranose, and AMP with AMPK protein, we could conclude that Asp151 was key residue for top 2 TCM candidates and the control bound with AMPK. By further analysis of distance of H-bonds between eugenyl beta-Dglucopyranoside and essential amino acids of AMPK, the compound formed H-bond with Asp151 only at early stage of MD. At late stage of MD, the compound formed H-bond with Gly22 instead. By further analysis of distance of H-bonds between 6-O-cinnamoyl-D-glucopyranose and essential amino acids of AMPK, the compound formed H-bond with Asp151 at all stages of MD. By further analysis of distance of H-bonds between AMP and essential amino acids of AMPK, the compound formed H-bond with Asp151 at all stages of MD. Interestingly, H36 of eugenyl beta-D-glucopyranoside formed H-bond with OD1 and OD2 of ASP151; H36 of AMP also formed H-bond with OD1 and OD2 of ASP151, but their figure patterns were quite different. However, H23 of 6-O-cinnamoyl-D-glucopyranose formed H-bond with OD1 and OD2 of ASP151 and the figure pattern was similar to that of AMP. We could confirm that all the candidates and the control formed stable complexes with AMPK. In addition, all the candidates and the control induced conformational change of AMPK due to changing position of H-bonds.

Conclusion
AMPK was a hot issue in the past decade because it plays an important role in sensing cellular energy or nutrient status. AMPK can regularize cell growth by the way of pseudostarvation leading to autophagy. Thus AMPK has become a great target for drug design of tumor inhibition. Activation of AMPK followed its conformational change. We tried to find potential compounds that could bind to AMPK by virtual screening of the world's largest TCM Database (http://tcm.cmu.edu.tw/) and had the ability to activate the function of AMPK. We selected eugenyl beta-D-glucopyranoside and 6-O-cinnamoyl-D-glucopyranose as candidates for further investigation. Through the methods of MD simulation consisting of RMSD, Gyrate, SAS, RMSF, total energy, DSSP, smallest distance matrices, occupancy and distance of H-bonds, and ligand pathway, we could conclude that the candidates and the control (AMP) formed stable complexes with AMPK. Asp151 was key residue for the candidates and the control bound with AMPK. These compounds also had the ability to induce AMPK conformational change. Thus eugenyl beta-D-glucopyranoside and 6-O-cinnamoyl-D-glucopyranose had the qualification to be AMPK agonist.