Lys169 of Human Glucokinase Is a Determinant for Glucose Phosphorylation: Implication for the Atomic Mechanism of Glucokinase Catalysis

Glucokinase (GK), a glucose sensor, maintains plasma glucose homeostasis via phosphorylation of glucose and is a potential therapeutic target for treating maturity-onset diabetes of the young (MODY) and persistent hyperinsulinemic hypoglycemia of infancy (PHHI). To characterize the catalytic mechanism of glucose phosphorylation by GK, we combined molecular modeling, molecular dynamics (MD) simulations, quantum mechanics/molecular mechanics (QM/MM) calculations, experimental mutagenesis and enzymatic kinetic analysis on both wild-type and mutated GK. Our three-dimensional (3D) model of the GK-Mg2+-ATP-glucose (GMAG) complex, is in agreement with a large number of mutagenesis data, and elucidates atomic information of the catalytic site in GK for glucose phosphorylation. A 10-ns MD simulation of the GMAG complex revealed that Lys169 plays a dominant role in glucose phosphorylation. This prediction was verified by experimental mutagenesis of GK (K169A) and enzymatic kinetic analyses of glucose phosphorylation. QM/MM calculations were further used to study the role of Lys169 in the catalytic mechanism of the glucose phosphorylation and we found that Lys169 enhances the binding of GK with both ATP and glucose by serving as a bridge between ATP and glucose. More importantly, Lys169 directly participates in the glucose phosphorylation as a general acid catalyst. Our findings provide mechanistic details of glucose phorphorylation catalyzed by GK, and are important for understanding the pathogenic mechanism of MODY.


Introduction
Glucokinase (GK) is a glycolic enzyme that catalyzes the phosphorylation of glucose to glucose-6-phosphate in the first step of glycolysis. Expressed in the liver, pancreas, brain and gut, GK plays a key role in maintaining glucose homeostasis through regulation of glucose-dependent insulin secretion in pancreatic cells, and glucose uptake and storage in the liver [1][2][3]. Clinical evidences have shown that GK might be an important therapeutic target for treating metabolic diseases such as maturity-onset diabetes of the young (MODY) and persistent hyperinsulinemic hypoglycemia of infancy (PHHI) [4][5][6][7].
Although GK belongs to the hexokinase family [8], GK is distinctive from other hexokinases because of its relatively low glucose affinity in the range of blood-glucose levels and its positive cooperative kinetics. These properties allow GK to rapidly respond to changes in glucose concentrations under physiological conditions, and thereby function as a glucose sensor [9]. Various studies have demonstrated that GK is a monomeric enzyme whose allosteric mechanism is different from those found in other extensively studied oligomeric allosteric enzymes [10][11][12][13][14]. Comparison of the crystal structures of GK in closed state (active state) and super-open state (inactive state) shows that GK may undergo a global conformational change between these two states, and such a conformational change may be responsible for the special allosteric characteristics of GK [15]. Recently, we used molecular dynamics simulation method to show that the global conformational transition pathway between the two states of GK includes three intermediate steps, identified by free energy scanning of snapshots throughout the pathway. The computational predictions were verified by mutagenesis and enzymatic kinetic analysis [16]. These studies gave the underlying principle of the enzymatic mechanism of GK and provided the explanations for the sigmoidal kinetic effect and the mnemonical mechanism for the cooperativity of GK with respect to glucose phosphorylation [17].
Recently, data for over 100 different mutations in the GK genes, potentially contributing to the development of an autosomal dominant form of type 2 diabetes, have been collected [18]. Among the data, some mutations in the MODY phenotype were found to be linked to alterations in both GK kinetics and regulation activities [19]. The catalytic mechanism of GK holds the most critical information for understanding the relation between the mutations in the MODY gene and GK regulation activity. This fundamental question remains unanswered, however, due to the lack of an accurate structural model for the complex of GK with ATP, glucose and Mg 2+ (GK-Mg 2+ -ATP-glucose complex). In particular, an important residue in the binding pocket, Lys169, can undergo mutation (K169N) in MODY [20] and its functional role remains unknown.
Here, we report a study of the catalytic mechanism of GK by combining molecular modeling, molecular dynamics (MD) simulations, quantum mechanics/molecular mechanics (QM/ MM) calculations, experimental mutagenesis and enzymatic kinetic analysis. The atomic structure of the glucose phosphorylation catalyzed by GK in aqueous solution was studied using protein substrate conformations obtained from the MD simulation. The simulation results are in agreement with the recent findings of mutagenesis experiments and related kinetic studies. MD simulations further revealed that Lys169 may play an essential role in both ligand binding and GK catalytic process. The computational prediction was verified by additional experimental mutagenesis and kinetic analysis. Based on these results, we propose an atomistic catalytic mechanism of GK for glucose phosphorylation using QM/MM calculations. Findings from this work provide a better understanding of the enzymatic mechanism of GK, and also a potential explanation of the pathogenic mechanism of MODY caused by the mutation in GK.

Results
The main goal of this study is to investigate the catalytic process of human GK by using computational methods in conjunction with enzymatic assay. Such a study, however, needs an accurate structural model for the whole catalytic environment of GK, i.e. the three-dimensional (3D) structure of the complex of GK with ATP, glucose and Mg 2+ (GK-Mg 2+ -ATP-glucose complex, designated as GMAG complex hereinafter). Accordingly, a 3D structural model of GMAG complex was constructed based on the crystal structures of GK and several other hexokinases by using molecular modeling. Following this, the dynamic conformation change of the catalytic site of GK was investigated by using MD simulations. Furthermore, the modeling and simulation results were validated using mutagenesis and enzymatic kinetic assays. Encouraged by the compatible results from theoretical prediction and experiments, we explored the GK catalytic mechanism for the glucose phosphorylation based on a stable conformation of GMAG complex derived from the MD simulations by using quantum mechanics/molecular mechanics (QM/MM) method.

3D Model of the GK-Mg 2+ -ATP-glucose (GMAG) Complex
The crystal structure of human GK had not been available until 2004 due to the high flexibility of the protein [21]. Moreover, it is easy for ATP to transfer its terminal phosphate to glucose in situ, making it difficult to gain accurate structural information of the catalytic environment of glucose inside GK, which plays a critical role in understanding the catalytic and regulation mechanisms of GK. Several models for GK catalytic environment have been published [22][23][24][25][26][27]. However, these models were built using homology modeling based on the X-ray crystal structure of human hexokinase type I [28,29]. Although human hexokinase type I is homologous to GK with sequence identity of ,55%, its three-dimensional (3D) structure is different, especially around the ATP binding pocket, as indicated by the structural superposition of the crystal structures of these two enzyme ( Figure S1). The reliability of these GK catalytic models is therefore questionable. On the basis of the crystal structures of GK, now available, we built up a structural model for the GMAG complex (see Materials and Method section for detail). This model is more accurate than the previous models, and has also been validated by further MD and QM/MM calculations, as well as mutagenesis experiments and enzymatic assays (see discussion below).
Structurally, GK is composed of a large domain and a small domain; the channel-shaped active site for phosphorylation is located in the deep cleft between these two domains, as shown in Figures 1A and 1B. The Mg 2+ ion and ATP (Mg 2+ -ATP) are predicted to bind to the left portion of the active site cleft, interacting with both domains [22][23][24]. The Mg 2+ -ATP binding site is formed by residues 78-83, 169, 225-228, 295-298, 331-336, 410-413, and 415-417 ( Figure 1C). Meanwhile, the glucose resides in a small pocket composed of residues 80-81, 151-154, 168-169, 204-207, 225-232, 254-259, 287 and 290 ( Figure 1D). Our 3D model of GMAG complex reveals that GK provides a favorable microenvironment for the phosphorylation of glucose. Several important hydrogen bonds between ligands (ATP, Mg 2+ and glucose) and GK are observed in our model (Table S1), which are in agreement with a large number of mutagenesis analysis data [20,[30][31][32][33][34]. The c-phosphate of ATP is close to the -O 6 H group of glucose with electrostatic interaction, indicating that the cphosphate is about to be transferred to glucose.
To study the stability of the reaction environment, a 10-ns molecular dynamics (MD) simulation was performed on the GMAG complex. The time evolutions of the weighted Root-Mean Square Deviations (wRMSD) for the atoms of GK, ATP and glucose from their initial positions (t = 0 ps) were monitored. The result indicates that all the reaction components are relatively stable except for the rotation of the -O 6 H group which causes an increase in wRMSD of glucose by 0.3 Å at 1 ns ( Figure S2). The H-bond network amongst GK, ATP and glucose also reflects the stability of the phosphorylation environment. These H-bonds were mostly maintained during the 10-ns MD simulation as indicated by their occupancies (Table S1).
Under the restriction of GK, ATP, Mg 2+ and glucose form a highly advantageous configuration for the biochemical reaction ( Figure 2 Figure 2 and Table S1). Most importantly, the -O 6 H group of glucose directly hydrogen bonds to the O c2 atom of ATP, producing the requisite complex for glucose phosphorylation. Residues detected in the naturally occurring mutations (e.g. K169N, T228M, E256K, S336L, and S411F) from MODY families [20,[30][31][32] are involved in the Hbond network for the binding of GK with glucose and ATP. Extensive work confirmed that S336L mutation decreased both the binding affinity of GK to ATP and the catalytic activity of the enzyme [33]. Several other site-directed mutants, including N204A, E256A and E290A, have been shown to alter GK's enzymatic kinetics by decreasing the binding affinity of glucose and lowering V max in the catalytic process [31,34]. These mutagenesis and enzymatic results demonstrate that our 3D model of GMAG complex is reasonable.

Importance of Lys169 in the Phosphorylation of Glucose and Mutation Validation
Among the conserved residues hydrogen bonding to ATP and glucose, Lys169 is of special interest because its cationic end, -N f H + 3 , bridges ATP and glucose together through H-bonds, placing the c-phosphorus (P c ) atom only 3.3 Å far way from the -O 6 H group of glucose ( Figure 2). This implies that Lys169 might play an important role in the phosphorylation of glucose. Additionally, it was reported that Lys169Asn (K169N) is one of the naturally occurring mutations in the GCK gene associated with familial mild fasting hyperglycemia [20]. The K169N mutant was also shown to experience a partial loss of glucose binding [30].
To address the biological function of Lys169, we constructed a 3D structural model for the GK K169A mutant (GK K169A ) in complex with ATP-Mg 2+ and glucose, and two additional MD simulations of the ATP-Mg 2+ -GK K169A and glucose-GK K169A complexes were conducted. Binding free energies of ATP and glucose with wild-type GK and GK K169A (Table S2) were then calculated using the MM-PBSA method encoded in AMBER (Version 8.0). Both ATP and glucose are able to bind tightly to the wild-type GK with calculated binding free energies of 218.6764.59 and 244.0663.94 kcal/mol, respectively. While neither ATP nor glucose can bind with GK K169A .
The MD simulations and binding free energy calculations indicate that Lys169 may dominate the binding of both ATP and glucose, and hence the K169A mutant possibly loses the catalytic function for the phosphorylation of glucose. To verify this conclusion, a bioassay on the K169A mutant was performed and result is shown in Table 1. The kinetic assay indicates that K169A is indeed an inactivating mutation; its enzymatic activity is completely abolished. Far-UV CD spectra and fluorescence emission spectra results demonstrated that the K169A mutant had well-defined secondary structures; thereby we can exclude the possibility that this mutation may have caused GK to undergo misfolding (Figures S3, S4, S5). Thus, the MD prediction, regarding the importance of Lys169 to the binding of ATP and glucose with GK and to the catalytic activity of GK as well was validated.

Lys169 Acts as a General Acid Catalyst
To further figure out the role of Lys169 in enzymatic catalysis, the reaction mechanism of glucose phosphorylation inside the active site of GK was investigated by using the QM/MM approach. The system for QM/MM simulation was constructed based on the snapshot at 2500 ps of the MD trajectory. The QM region was composed of Mg 2+ , water molecules coordinating with Mg 2+ , important groups of ATP, glucose, Lys169 and Asp205; the remainder of GK was included in the MM region ( Figure S6). The partitioning scheme for QM and MM regions is described in the Materials and Methods section and Figure S6. We designate this structure as reagent system. For the reagent system, the initial coordination configuration of Mg 2+ with waters, ATP and glucose and the hydrogen bonding pattern amongst ATP, glucose, Lys169 and Asp205 have been described above ( Figure 2). ONIOM, a QM/MM method encoded in Gaussion03 [35], was used for all the QM/MM calculations (Figure 3 and 4).
The QM/MM optimized geometry of the reagent system showed that Mg 2+ octahedrally coordinates with the O b1 and O c1  atoms of ATP, the O d2 atom of Asp205 and three oxygen atoms from three water molecules ( Figure 3). The major difference of the optimized structure for the reaction center from the initial structure obtained by MD simulation is that one water molecule substituted the O c3 atom to coordinate with the Mg 2+ ion. QM optimization resulted in important structural changes viz. one proton of the -N f H + 3 group of Lys169 transferred to the O c3 atom, and the resulting -N f H 2 group formed a H-bond with the new water molecule coordinating with the Mg 2+ ion; on the other hand, the -O 6 H group of glucose adjusted its direction to point toward the O d1 atom of Asp205, forming two strong H-bonds with the O d1 and O d2 atoms of Asp205. This structural reorganization made the -O 6 H group more propitious to attack the P c atom of ATP, implying that Lys169 might act as an acid catalyst and Asp205 as a base catalyst in the phosphorylation of glucose. Thus, we propose a mechanism for the glucose phosphorylation catalyzed by GK as Scheme 1 ( Figure 5). Along this reaction path, the energies of the reagent (R), transition state (TS), and immediate product (P) were determined by two-dimensional QM/ MM potential energy surface by defining the distances of R(O b3 _P c ) and R(P c _O 6 ) as the reaction coordinates ( Figure 3). In the optimized reagent, R(O b3 _P c ) = 1.69 Å and R(P c _O 6 ) = 3.39 Å ; and in the optimized immediate product, the O b3 _P c bond is broken and the P c _O 6 is formed, R(P c _O 6 ) = 1.71 Å . The calculated potential energy barrier of Scheme 1 is DE ? = 18.3 kcal/mol.
The structure of the transition state (TS) in Scheme 1 was determined by adiabatic mapping at the QM/MM level (Figure 3). In the TS, R(O b3 _P c ) = 2.1 Å and R(P c _O 6 ) = 2.0 Å , and the cphosphate (-P c O 3 ) is turned to be in a plane. This structure clearly illustrates that the P c _O 6 bond is partially formed and the O b3 _P c bond is partially broken. The overall reaction is calculated to be exothermic by DE = 222.1 kcal/mol.
To further address the catalytic role of Lys169, we investigated the reaction mechanism of glucose phosphorylation catalyzed by the GK K169A mutant employing the QM/MM method. Compared with the optimized structure of GK reagent system, the major difference for the ATP-GK K169A -glucose reaction center is that the -O 6 H group of glucose doesn't form H-bond with either O d1 or O d2 atom of Asp205, but forms a strong H-bond with the O c3 atom of ATP. This optimized structure implies an alternative path for the glucose phosphorylation as Scheme 2 ( Figure 6): the proton of the -O 6 H group transfers to the O c3 atom of ATP while the O 6 atom attacks the P c atom of ATP, and at the same time the O b3 _P c bond is broken. According to this reaction path, we obtained the structures and energies of reagent (R), transition state (TS), and immediate product (P) by calculating the twodimensional QM/MM potential energy surface; the distances of R(O b3 _P c ) and R(P c _O 6 ) were defined as the reaction coordinates. The result is shown in Figure 4. For the TS, R(O b3 _P c ) = 2.1 Å and R(P c _O 6 ) = 1.9 Å . Remarkably, the calculated potential energy barrier of this reaction path is DE ? = 32.1 kcal/mol, which is ,14 kcal/mol higher than that of the reaction catalyzed by the wild-type GK. This result indicates that the phosphorylation of glucose could not occur without the assistance of Lys169.

Discussion
As mentioned above, GK is a key enzyme that phosphorylates glucose and triggers glucose utilization and metabolism. As an attractive drug target for discovering anti-diabetes drug, the importance of this enzyme has been appreciated. However, the unique catalytic mechanism of GK still remains unclear [17,18]. Several mutants of GK with the MODY phenotype have been identified recently; in particular, it was found that K169N is a naturally occurring mutation (K169N) in MODY [20]. In this study, we have investigated the phosphorylation mechanism of glucose catalyzed by GK, and revealed that Lys169 is a crucial residue for glucose phosphorylation.
By using molecular modeling and simulation methods, we constructed a 3D structural model for the complex of GK with ATP, glucose and Mg 2+ (GMAG complex). This structural model is more reliable than previously published models [22][23][24][25][26][27]36] because it was constructed based on the crystal structures of human GK. In particular, this model provides a more realistic reaction environment for the phosphorylation of glucose. The validity of the structural model of GMAG was confirmed by the agreement between the available experimental mutagenesis and enzymatic data for GK and the important interactions of GK with ATP, glucose and Mg 2+ identified through computation [20,[30][31][32][33][34].
Both ATP and glucose interact with the residues lining the binding pocket of GK through H-bonds (Figure 2). Among these residues, Lys169 plays a significant role, because it forms H-bonds with both ATP and glucose. This indicates that Lys169 contributes to the binding of GK with either ATP or glucose or both. Computational simulations demonstrated that both ATP and glucose are able to bind to the wild-type GK, and neither ATP nor glucose can bind to the GK K169A mutant (Table S2). Indeed, the experimental mutagenesis and enzymatic kinetic assay validated the computational prediction, i.e. the mutant of K169A could not bind to ATP and glucose and its enzymatic activity is completely lost ( Table 1).
The QM/MM calculations on the mechanisms of glucose phosphorylation catalyzed by both GK and GK K169A mutant further revealed the role of Lys169 in catalysis. During the phosphorylation process of glucose, in addition to enhancing the binding of ATP and glucose with GK, Lys169 also acts as a general acid catalyst, providing a proton to protonate the O c3 atom of ATP; at the same time Asp205 acts as a general base catalyst, extracting the proton from the -O 6 H group of glucose (Scheme 1). On the other hand, ATP and glucose bind to the GK K169A mutant with a different model than with the wild-type GK (Figures 3 and 4), whereby the phosphorylation of glucose adopts a different mechanism (Scheme 2). The QM/MM calculated activation energy of Scheme 2 is ,14 kcal/mol higher than that of Scheme 1. This suggests that the GK K169A mutant has significantly reduced catalytic activity than its wild-type form. The calculated result is in good agreement with our mutagenesis data (Table 1). Moreover, the QM/MM result also indicates that Asp205 of GK loses its function as a general base catalyst in the glucose phosphorylation when Lys169 is absent (Figure 4).
The agreement between computation and experimental results validates the essential role of Lys169 in the phosphorylation of glucose. In summary, Lys169 plays at least three functional roles in the metabolism of glucose: (1) it enhances the binding of GK to both ATP and glucose; (2) it bridges ATP and glucose together; and (3) it acts as a general acid catalyst and directly participates in glucose phosphorylation. These results are important for understanding the catalytic mechanism of GK and the cause of the pathogenic mechanism of MODY due to the GK mutation [17,18]. Obviously, asparagine cannot act as a general acid catalyst because its side chain is incapable of providing a proton to ATP, as shown in Figure S7. Our study thus provides an explanation of why the catalytic activity of the mutant GK (K169N) in the MODY is lower than that of the wild-type GK.

Structural Model
On the basis of crystal structures GK (the active human GK with glucose, PDB entry 1V4S) and two hexokinases (human brain hexokinase complexed with glucose and phosphate, PDB entry 1HKC, and monomeric hexokinase I complexed with ADP, PDB entry 1DGK) as templates, the whole structural model of GMAG complex (GK-Mg 2+ -ATP-glucose complex) was constructed by using a heuristic approach. The crystal structure of GK-glucose complex (1V4S) was used as a basal structure. The ATP conformation was generated according to the structure of ADP in the crystal structure of 1DGK, and the coordination configuration of Mg 2+ with ATP (Mg 2+ -ATP) was constructed based on the coordination state of Mg 2+ in the crystal structure of 1HKC. Then the Mg 2+ -ATP was ''grafted'' to the binding pocket of GK, and the primary position of Mg 2+ -ATP was decided through structural superimposition for the crystal structures of the three enzymes using the Homology module of Insight2005 (Accelrys, San Diego, CA). For the initial structure of the GMAG complex, residues within a radius of 6.5 Å around the binding site were optimized by energy minimization using the AMBER force Lys169 Plays a Key Role in GK PLoS ONE | www.plosone.org field implemented in the Sybyl6.8 (Tripos, St. Louis, MO) with following parameters: a distance-dependent dielectric function, non-bonded cutoff of 8 Å , Amber charges for the protein, Gastieger-Hückel charges for ATP and glucose. The structure was minimized by simplex method first, followed by Powell method to an energy gradient,0.05 kcal/(mol?Å ). The structural model of G K169A MAG and G K169N MAG complex were simply constructed by replacing Lys169 with Ala and Asn respectively, followed by minimization using the same protocol as in the structural optimization for the model of GMAG.

Molecular Dynamics Simulation and Binding Free Energy Calculation
MD simulations were performed on the GMAG and G K169A -MAG complexes. Before simulations, each complex was put into a suitably sized box, of which the minimal distance from the protein to the box wall was 15 Å . Then the box was solvated with the SPC (Simple Point Charge) water model [37]. Each complex/ water system was submitted to energy minimization. Afterwards, counterions were added to the system to provide a neutral simulation system. The whole system was subsequently minimized again. The charges of the atoms of ATP and glucose were calculated by using the RESP method [38] encoded in the AMBER suite of programs [39] at the level of RHF/6-31G*. Covalent and non-bonded parameters for the ligand atoms were assigned, by analogy or through interpolation, from those already present in the AMBER force field [40].
MD simulations were carried out using the AMBER package (Version 8.0) with constant temperature, constant pressure (NPT), and periodic boundary conditions. The AMBER Parm99 force field [40] was applied for the proteins. The Particle Mesh Ewald method [41] was used to calculate the long-range electrostatics interactions. The non-bonded cutoff was set to 12.0 Å , and the non-bonded pairs were updated every 25 steps. The SHAKE method [42] was applied to constrain all covalent bonds involving hydrogen atoms. Each simulation was coupled to a 300 K thermal bath at 1.0 atm of pressure (1 atm = 101.3 kPa) by applying the algorithm of Berendsen et al. [43]. The temperature and pressure coupling parameters were set as 1 ps. An integration step of 2 fs was set up for the MD simulations.
On the basis of the equilibrated dynamic trajectory, the binding free energy of each system (GMAG, GK K169A -Glucose and GK K169A -Mg 2+ -ATP) was calculated by using the MM-PBSA method encoded in the AMBER 8.0 program. Coordinates from the trajectory were used every 20 ps, and the MM-PBSA calculation was performed on each of them using the AMBER 8.0 program. For each snapshot collected during the simulations, binding free energy (DG binding ) was calculated using Eq. (1): where DG complex , DG ligand and DG receptor are the free energies of the complex and two components of interaction. Each free energy term in Eq. (1) in AMBER 8.0 was calculated with the absolute free energy in gas phase (E gas ), the solvation free energy (DG solvation ) and the entropy term (TDS) using Eq. (2): E gas is the sum of the internal strain energy (E int ), van der Waals energy (E vdw ) and electrostatic energy (E electrostatic ) (Eq. (3)). E int is the energy associated with vibrations of covalent bonds and bond angles, rotation of single bond torsional angles (Eq. (4)).
E gas~Eint zE vdw zE electrostatic ð3Þ The solvation free energy, DG solvation , is approximated as the sum of the polar contribution (G PB ) and nonpolar contribution (G nonpolar ) using a continuum representation of the solvent: The polar contribution (G PB ) to the solvation energy was calculated using PB model in sander module. The nonpolar contributions (G nonpolar ) were estimated using a simple equation: G nonpolar = c6SASA+b (kcal/mol). SASA is the solvent-accessible surface area estimated using the MSMS algorithm with probe radius of 1.4 Å . The surface tension proportionality constant c and the free energy of nonpolar solvation for a point solute b were set to 0.00542 kcal/(mol Å 2 ) and 0.092 kcal/mol, respectively.
The entropy calculation is extremely time-consuming for large systems. In addition, the main aim of calculating the binding free energy was to address the influence of mutations to the binding affinity. Site-directed mutation may not result in dramatic conformational changes for the protein and it was assumed that entropy contributes little to the relative binding free energy for the binding of mutant in comparison with the wild-type protein. Therefore, in this study, DG binding (without term of -TDS) was estimated to address the mutation effect to of the binding free energy.

QM/MM Calculation
QM/MM calculations were performed by using a two-layered ONIOM method encoded in the Gaussian03 program [35]. The ONIOM method is a hybrid computational method developed by Morokuma and coworkers that allows different levels of theory to be applied to different parts of a molecular system [44][45][46][47][48][49][50][51][52]. In the two-layered ONIOM method, the molecular system under study is divided into an inner layer and an outer layer. The inner layer consists of the most critical elements of the system, and the rest of the system comprises the outer layer. In the terminology of Morokuma and co-workers, the full system is called ''real'' and is treated with a low level of theory. The inner layer is termed ''model'' and is treated with both the low level of theory and a high level of theory. The total ONIOM energy E ONIOM is given as following: where E(high, model) is the energy of the inner layer (plus the link atoms) at the high level of theory, E(low, real) is the energy of the entire system at the low level of ONIOM theory, and E(low, model) is the energy of the model system at the low level of theory. Thus, the ONIOM method allows one to perform a high-level calculation on just a small, critical part of the molecular system and incorporate the effects of the surrounding elements at a lower level of theory to yield a consistent energy expression on the full system. The minimized structure optimized using the AMBER Parm99 force field was further optimized at the ONIOM (B3LYP/6-31G*: Amber) level. The quantum mechanics (QM) region includes the b, c-phosphate group of fully unprotonated ATP, four water molecules including two coordinated water molecules, Mg 2+ ion, the hydroxymethyl group (-CH 2 -O 6 H) of glucose, the methyleneam-monium group (-CH 2 -N f H + 3 ) of Lys169 and the anionic carboxymethyl group (-CH 2 -COO 2 ) of Asp205 for a total of 40 atoms ( Figure S6). The QM region was calculated by using the density functional theory with the B3LYP exchange-correlation functional and 6-31G* basis set. The remainder of the system (MM region) was treated by using the AMBER Parm99 force field. A total of 7,147 atoms were included for the QM/MM calculations by using the ONIOM module as implemented in Gaussian 03. The electrostatic interactions between the QM and MM regions were calculated by using the electronic embedding method, which treats the polarization of the QM region by the MM region with scaled partial atomic charges of MM atoms, and the response of QM region with Merz-Singh-Kollman scheme for charge fitting so as to produce the changing partial charges of the QM atoms. The charge for the QM region was 21, and the charge for the MM was 1, therefore the total system remained neutral.

Expression, Purification and Kinetic Analysis on GK and Its K169A Mutant
Escherichia coli (E. coli) host strain M15 was purchased from Qiagen (Valencia, CA). All chemicals were of reagent grade or ultra-pure quality, and commercially available. Human liver glucokinase isoform 2 (hLGK2) DNA sequence was prepared by total gene synthesis method (Sangon, Shanghai, China).
Expression and purification of the recombinant protein was performed according to His-Bind Kits (Novagen, San Diego, CA). The recombinant clones pQE-30-hLGK2 and pQE-30-hLGK2 K169A were separately transformed into E. coli strain M15 that grew in LB media supplemented with 100 mg/mL ampicillin and 50 mg/mL kanamycin at 37 uC. Expressions of pQE-30-hLGK2 was induced by 0.3 mM IPTG (Isopropyl-b-D-thiogalactopyranoside) when the OD 600 reached 0.6 at 25 uC for an additional 16-20 h. The hLGK2 cells were harvested by centrifugation and suspended in Buffer A (20 mM Tris-HCl, pH 8.0, 500 mM NaCl, 10 mM imidazole), stored at 280 uC after centrifugation until use. Packed cells were suspended in Buffer A plus 1 mM PMSF (Phenylmethanesulfonylfluoride), after sonication treatment on ice, the mixture was centrifuged to yield a clear supernatant, which was loaded onto a column with Ni-NTA resin (Qiagen) preequilibrated in Buffer A (20 mM Tris-HCl, pH 8.0, 500 mM NaCl, 10 mM imidazole). The column was washed with 8-10 column volumes of Buffer B (20 mM Tris-HCl, pH 8.0, 500 mM NaCl, 20 mM imidazole) and eluted with Buffer C (20 mM Tris-HCl, pH 8.0, 500 mM NaCl, 150 mM imidazole), then the soluble hLGK2 fractions were pooled separately and dialyzed against Buffer D (25 mM HEPES, pH 7.1, 25 mM KCl, 2 mM MgCl 2 , 1 mM DTT) to remove imidazole. The hLGK2 protein was thus concentrated by ultrafiltration with Amicon centrifugal filter device to appropriate concentration. All procedures including purification, dialysis and concentration were performed at 4 uC. Protein concentration was determined by Bradford assay using bovine serum albumin (BSA) as standard. The expressions, purifications and the subsequent procedures of hLGK2 mutant K169A protein were the same as that of hLGK2 protein.
The recombinant proteins of hLGK2 and hLGK2 mutant K169A was subjected to kinetic analysis (V max , K m-ATP , S 0.5 and nH) according to the published procedure reported by Grimsby et al. [53]. CD spectra of hLGK2 WT and hLGK2 mutant proteins CD spectra of hLGK2 WT and K169A mutant proteins were recorded on a Jasco-810 spectropolarimeter equipped with a thermal controller. The protein sample was prepared in 20 mM NaH 2 PO 4 , pH 7.1, at 25uC with protein concentration of 5 mM. Far-UV CD spectra from 190 to 250 nm were collected with 1 nm bandwidth using a 0.1 cm path length cuvette and normalized by subtracting the base line recorded for the buffer under identical conditions. Each measurement was repeated three times, and the final result was the average of the three independent scans.

Fluorescence Spectra of Wild-type hLGK2 and hLGK2 Mutant Proteins
Fluorescence experiments of hLGK2 WT and mutant hLGK2 proteins K169A were carried out on a Hitachi F-2500 fluorescence spectrophotometer. The samples were prepared in 20 mM NaH 2 PO 4 solution at pH 7.1 with protein concentration at 5 mM. The fluorescence emission spectra from 300 to 400 were collected after excitation at 280 nm. Fluorescence emission spectra of proteins were determined using a 1 cm path quartz cuvette at 25uC. The spectral band width was 5 nm for excitation and 10 nm for emission. The maximal emission wavelength was calculated by the affiliated software of the F-2500 spectrophotometer.