A Description of Enzymatic Catalysis in N-Acetylhexosamine 1-Kinase

41 mins ago - Herein, by using state-of-the art ab initio QM/MM MD and MM MD simulations, one-dimensional and two-dimensional free energy profiles to ...
0 downloads 7 Views 4MB Size
Subscriber access provided by UNIV OF CAMBRIDGE

A Description of Enzymatic Catalysis in N-Acetylhexosamine 1Kinase: Concerted Mechanism of Two-Magnesium-Assistant GlcNAc Phosphorylation, Flexibility Behavior of Lid Motif upon Substrate Recognition, and Water-Assistant GlcNAc-1-P Release Yuan Zhao, Nai She, Yiming Ma, Chaojie Wang, and Zexing Cao ACS Catal., Just Accepted Manuscript • DOI: 10.1021/acscatal.8b00006 • Publication Date (Web): 04 Apr 2018 Downloaded from http://pubs.acs.org on April 4, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

A Description of Enzymatic Catalysis in N-Acetylhexosamine 1Kinase: Concerted Mechanism of Two-Magnesium-Assistant GlcNAc Phosphorylation, Flexibility Behavior of Lid Motif upon Substrate Recognition, and Water-Assistant GlcNAc-1-P Release Yuan Zhao†, Nai She†, Yiming Ma§, Chaojie Wang*†, and Zexing Cao*‡ †

The Key Laboratory of Natural Medicine and Immuno-Engineering, Henan University, Kaifeng 475004, People’s Republic of China ‡

State Key Laboratory of Physical Chemistry of Solid Surfaces and Fujian Provincial Key Laboratory of Theoretical and Computational Chemistry, College of Chemistry and Chemical Engineering, Xiamen University, Xiamen 360015, People’s Republic of China §

College of Chemistry and Chemical Engineering, Henan University, Kaifeng 475004, People’s Republic of China

ABSTRACT: The N-acetylhexosamine 1-kinase (NahK) is the typical example of anomeric kinases acting on gluco-type substrate, which catalyzes the phosphorylation of GlcNAc or GalNAc at anomeric C1 position with ATP, playing a crucial role in bifidobacteria metabolic pathway and biosynthesis of sugar 1-phosphates and oligosaccharides. Herein, by using state-of-the art ab initio QM/MM MD and MM MD simulations, one-dimensional and two-dimensional free energy profiles to descript catalytic process have been explored. A concerted mechanism has been recognized for the delivery of phosphate group and proton to product ADP and GlcNAc-1-P with the free energy barrier of ~ 7.0 kcal/mol. The chemical reaction is assisted by two magnesium ions with a standing six coordination structure during the enzymatic process. The deficiency of Mg2+ bridging β-γ phosphates of ATP may alter the catalytic mechanism and brings higher free energy barrier. The protonation of GlcNAc-1-P is benefit for its getting rid of the magnesium ion binding. The water-assistant GlcNAc1-P cleavage from binding of Mg2+ needs ~ 9.4 kcal/mol at least, which means that the Pi-Pi bond breaking in chemical reaction step is probably not the rate-limiting step in the whole enzymatic process. A strongly exothermic phenomenon and an open-closed structural change of lid motif have been observed upon the GlcNAc binding. The exothermic trend strongly depends on the quantity and quality of hydrogen bond network around the ligand. KEYWORDS: N-acetylhexosamine 1-kinase, magnesium ion, catalytic mechanism, QM/MM MD, substrate binding

1.

INTRODUCTION

The bifidobacteria inhabited in the lower intestine of humans are considered as the probiotics agents, since it is benefit to human metabolism and immune systems, which play a crucial role in prevention and treatment of gastrointestinal disorders for humans and animals1-2. The digestible sugars are degraded and absorbed by the living beings during digestion, which are scarce in the lower intestines. Thus the bifidobacteria gain energy from indigestible oligosaccharides and glycoconjugates by using a variety of glycosidases, transporters, and metabolic enzymes3-4. A distinct lacto-N-biose (LNB)/galacto-N-biose (GNB) pathway has been identified in bifidobacteria that utilize LNB and GNB from human milk oligosaccharides (HMO) and glycoconjugates, which play an important role in colonization in the infant gastrointestinal tract. The LNB/GNB pathway is constituted by a specific ABC

transporter and four intracellular enzymes: galacto-Nbiose/lacto-N-biose phosphorylase (GLNBP), Nacetylhexosamine 1-kinase(NahK), UDP-glucose hexose 1phosphate uridylyltransferase (GalTBL), and UDP-glucose 4-epimerase (GalEBL), which is an energy saving variation of the Leloir pathway3, 5. Herein we mainly focus on NahK, which is the first example of anomeric kinases acting on a gluco-type substrate5. It can catalyze the phosphorylation of sugars Nacetyl-D-glucosamine (GlcNAc) or N-acetyl-Dgalactosamine (GalNAc) at the anomeric C1 position with ATP, as we can see in Figure 1a. The GlcNAc and GalNAc are the common primary metabolites in all living beings, and the product GlcNAc-1-phosphate (GlcNAc-1-P) or GalNAc-1-phosphate (GalNAc-1-P) will be further metabolized by the other enzymes to enter an amino-sugar metabolic pathway, which has important significance for

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 18

(a)

(b)

(c)

Figure 1. (a) NahK catalyzes the phosphorylation of GlcNAc or GalNAc at the anomeric C1 position with ATP. (b) The crystal 2+ 2+ structures of NahK•GlcNAc-1-P•ADP contained one Mg (PDB ID: 4OCP) and NahK•ATP contained two Mg (PDB ID: 4WH3). (c) The probable catalytic mechanism of phosphorylation of GlcNAc and GalNAc by NahK in experiment.

their host survival2, 5. It is worth mentioning that NahK has been used in enzymatic synthesis of oligosaccharides and sugar 1-phosphates for its broad substrate specificity towards structurally modified corresponding sugar analogue6-7; moreover, the cost of synthesis is reduced dramatically6-8. Therefore, clarification of the enzymatic catalytic mechanism is becoming extremely crucial in the understanding physiological pathology process and assist of various sugar’s biosynthesis. Experimentally, the crystal structures of NahK have been reported9-10. Wang et al.9 obtained seven threedimensional structures of NahK in complex with GlcNAc, GalNAc, GlcNAc-1-P, GlcNAc•AMPPNP, and GlcNAc-1P•ADP from both bifidobacterium longum and infantis at resolutions of 1.5 ~ 2.2 Å. Moreover, one Mg2+ has been

also identified. They found that NahK is a monomer in solution that shows a crescent-like structure with two domains to form a deep cleft (see Figure 1b). The crystals are grown in the presence of N-acetylhexosamine sugars and nucleotide ligands, and thus all the complex structures present a closed state. Besides that, they put forward a probable catalytic mechanism of phosphorylation of GlcNAc and GalNAc by NahK, as shown in Figure 1c. First of all, the substrate is deprotonated by Asp208, and then accepts the phosphate group from ATP. The Mg2+ is coordinated with α and β phosphates of ADP, and finally, the product sugar 1-phosphates and ADP are released. Herein, they conjectured the GlcNAc-1-P/GalNAc-1-P probably leaves prior to ADP on basis of crystal structures. The study provides the important information for understand-

ACS Paragon Plus Environment

Page 3 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

ing the NahK, while the enzyme is incorrectly named as “N-acetylhexosamine 1-phosphate kinases”, which should be corrected in the future related studies. Sato et al.10 determined an apo state of NahK and two complex crystals of NahK•ADP and NahK•ATP; moreover, all of them show an open state (see Figure 1b). They speculated that an induced fit structural change defined by two rigid domains. Most noteworthy is that the two magnesium ion has been discovered, bridging the α-γ and β-γ phosphates, which are essential for enzymatic catalysis by mutational analysis. They also mentioned that the action of the two magnesium ions is similar with protein kinases and aminoglycoside phosphotransferases, but distinct from the structures of other anomeric sugar kinases. Both experimental studies9-10 have speculated a molecular basis for the substrate recognition and catalytic mechanism for NahK. To sum up, the whole enzymatic process includes entrance of GlcNAc/GalNAc and ATP into the active site of protein, chemical reaction, and release of GlcNAc-1-P/GalNAc-1-P and ADP. However, mechanistic details for such processes are still lack theoretically. For example, the significant factor of substrate recognition and the mechanism of chemical reaction are not clear; the role of each stage in the whole process is not identified; moreover, the secondary Mg2+ in the active site structural similar protein always has more complex kinetic characteristic in the catalysis11, and what role does the secondary Mg2+ play in the NahK system is still open. Accordingly, we try to build a picture of the catalytic process, and thus the following crucial issues will be dealt with in the present study: (i) the dynamic and thermodynamic characterization of the GlcNAc transportation and phosphorylation; (ii) the possible channels for reactant delivery; (iii) the conformational changes of different domains in substrate recognition; (iv) the role of the two magnesium in the catalytic reaction; (v) the cleavage of GlcNAc-1-P from binding of Mg2+. To clarify these issues, extensive classical molecular dynamics (MD) and combined quantum mechanics and molecular mechanics (QM/MM) MD simulations have been employed. A thorough understanding of substrate binding and enzymatic mechanisms of NahK can be expected to be useful for further experimental and theoretical researches and enzymatic synthesis of various sugar molecules. 2.

COMPUTATIONAL METHODS

In the present work, we mainly focus on the NahK in complex with ATP and GlcNAc. Two X-ray crystal configurations are selected to build the initial simulation model: NahK•ATP complex in open state (PDB ID: 4WH3)10 and NahK•GlcNAc•AMPPNP complex in closed state (PDB ID: 4OCK)9. Compared with the crystal of 4WH3, the structure of 4OCK at high resolution has more complete amino acid sequence, and moreover the location of substrates is more clearly. However, 4WH3 has two magnesium ions that are crucial in catalytic process, and its open state is benefit for understanding the conformational changes of protein in substrate recognition and catalysis. Hence the

two initial models for the molecular dynamics simulations are constructed by employing the protein structures of both crystals (the protein of 4OCK is selected as Model A, and the protein of 4WH3 is chosen as Model B). Herein, the two magnesium ions of 4WH3 as well as the location of AMPPNP and GlcNAc for 4OCK are applied, and AMPPNP is changed into ATP. The overlapping scheme and preliminary structure adjustment have been performed. The protonation states of ionizable residues are determined by employing the program H++12-15 as well as their surrounding environment. The whole system is solvated into a cubic TIP3P16 water box of 89 × 85 × 76 Å with a 10 Å buffer distance on each side. Several Na+ are added to neutralize the model system. The ATP and GlcNAc are treated with Amber GAFF force field17, and the partial atomic charge is handled with the restrained electrostatic potential (RESP) charge18 by using HF/6-31G(d) calculation with Gaussian 09 package19. The force field parameters of protein are generated from Amber99SB force field20-21. The Lennard-Jones parameters of Mg2+ are disposed with the compromise parameter set from Li et al.22. The initial topology and coordinates are obtained by the tleap module in Amber 12 software23. For each model, multistep molecular dynamics simulations are carried out. Firstly, the entire system is minimized by three steps to remove poor interatomic contacts: (i) all the water molecules are optimized with constraining the enzyme, substrate, and two magnesium ions; (ii) the enzyme is optimized while keeping the substrate and magnesium ions restrained; (iii) the whole system is relaxed without any constrain force. Secondly, the system is heated up from 0 K to 300 K gradually under NPT ensemble for 100 ps, and then another 100 ps NPT MD simulations are carried out to relax the density to about 1.0 g/cm3. Finally, 20 ns NVT MD simulations are performed by employing the periodic boundary condition with 12 Å cutoff distances for electrostatic and van der Waals interactions calculations. Long-range electrostatic interactions are handled with particle mesh Ewald (PME) method24-25. The Langevin method has been applied to control the temperature at 300 K. The bonds with hydrogen atoms are constrained by using the SHAKE algorithm26. It is worth mentioning that when two models reach the stable state, they are almost overlapping and more likely to show a closed state. The two magnesium ions are only binding with ATP, and the GlcNAc has broken away from ATP. The distance between ATP and GlcNAc is suitable for GlcNAc transport study, while is not appropriate for the chemical reaction mechanism calculations. Thus another model of enzyme complex with ADP and GlcNAc-1P is constructed (named as Model C), which considers the effect of the two magnesium ions on both two substrates. The new system is also calculated according to the above simulation scheme. Another noteworthy is that an open and closed conformational changes will occur upon GlcNAc binding according to the experimental study, and hence the NahK complex with ATP is also calculated on

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2+

Page 4 of 18

2+

Figure 2. The QM/MM models of NahK•GlcNAc•ATP with two Mg (Model E) and one Mg (Model F).

basis of the above simulation scheme (named as Model D), which will be compared with the complex of NahK•ATP•GlcNAc. The possible channels of GlcNAc accessing into NahK are explored by combined random acceleration molecular dynamics (RAMD) and MD simulations (RAMD-MD)27-28. The snapshot from the stable trajectories of Model A is used as the initial structure. The protein and ligands are treated with the same force field parameters as molecular dynamics simulations. The random accelerations of 0.50, 0.40, 0.30, and 0.20 kcal Å-1 g-1 and a threshold distance of 0.005, 0.004, 0.003, and 0.002 Å are added to the center of mass of GlcNAc. If the ligand moves away from the threshold parameter within a certain period of time, the direction is maintained; otherwise a new random direction will be selected. When the ligand flees away from the initial location, the classic MD simulations will be started to relax the sampling, which can be avoided the wrong direction caused by the higher random force. Based on the results of RAMD MD simulations, the most favorable channel will be analyzed. The thermodynamics and dynamics properties as well as the protein conformational changes are identified by using MM MD simulations combined with umbrella sampling technology29. On basis of the appropriate biasing harmonic potential, 30 ns NVT MD simulations are carried out for each window along the defined reaction coordinate. The free energy profile is mapped by employing the weighted histogram analysis method (WHAM) technology30-31 with the probability distribution of each window. With the final snapshot from the stable trajectories of Model C, the QM/MM model is prepared by removing all water molecules beyond 20 Å from the second magnesium ion (named as Model E), which is shown in Figure 2. The QM subsystem includes the Asp208, Lys210, Asn209,

and Asp228, GlcNAc-1-P, ADP, two magnesium ions, and the coordinated water molecules, which are described by using the B3LYP functional32-33 and the 6-31G(d) basis set. This level has been employed successfully in other enzymes34-37. The other atoms are classified as the MM subsystem. The protein is calculated by employing the Amber99SB force filed20-21, and water molecules are treated with TIP3P model16. The QM/MM boundary is depicted by the improved pseudobond approach38-41. The spherical boundary condition is applied, and all the atoms are free to move. The cutoff radii of 18 and 12 Å are utilized for electrostatic and van der Waals interactions among MM atoms, respectively. With regard to the QM/MM system, the optimization is performed, and then the reaction coordinate driving (RCD) method42 is used to obtain the minimum reaction energy path (MEP); meanwhile, the reactant point is identified. In order to map the forward reaction picture, the reactant state is relaxed by using QM/MM MD simulations (see Figure 2), and then the RCD method is applied to acquire the new MEP. For all of the electronic structures along the reaction path, the MM region is equilibrated for 500 ps by MM MD simulations with QM subsystem frozen. The snapshots are serves as the initial structures for subsequent QM/MM MD simulations with umbrella sampling29 technology. Each window is simulated for 20 ps with 1 fs time step that applied the Beeman algorithm43. Herein, the Berendsen thermostat44 is used to control the system at 300 K. The probability distribution along the defined reaction coordinate is obtained, and last 4 ps trajectories for each window are collected to determine the free energy profile with the WHAM technology30-31. All the ab initio QM/MM simulations are carried out by using Q-Chem 4.045 and Tinker 4.246 packages.

ACS Paragon Plus Environment

Page 5 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

In order to investigate the role of the secondary Mg2+, the 20 ps QM/MM MD simulations are performed to relax the system based on Model E by removing the secondary Mg2+ that bridging the α-γ and β-γ phosphates of ATP. We found that the coordination number of Mg2+ is changed from six to four after simulations, which named as Model F. The QM subsystem includes Asp208, Lys210, Asn213, Asp228, ADP, GlcNAc-1-P, one magnesium ion, and coordinated water molecules, which are calculated at the B3LYP/6-31G(d) level. The other atoms are classified as MM subsystem, which are treated with the Amber99SB force field. The calculation process is same as Model E, which is carried out to obtain structures of key states and the thermodynamic characterization of the reaction. 3. RESULTS AND DISCUSSIONS 3.1. Comparison of Models A-D For each model, on basis of the root-mean-square deviation (RMSD) results of all the backbone atoms in Figure S1, we can see that after 8 ns MM MD simulations, the RMSD values are fluctuated steadily around 1.42, 1.76, 1.07, and 1.28 Å for Model A, B, C, and D, respectively, and their SDs are 0.12, 0.14, 0.08, and 0.11 Å, which means that the trajectories reach equilibrium. The snapshots of stable state are abstracted for the protein conformational analysis (see Table S1). With regard to Model A, as is shown in Figure 3, its initial state and stable state after MD simulations almost overlap, which show a closed state; for Model B, an open state is observed for its original crystal, while after MD simulations, the “lid motif” that refers to loop 26 ~ 32 and loop 314 ~ 321 displays a tendency of closed state. It means that the NahK is inclined to present a closed state, when both ATP and GlcNAc are binding in the active site. For Model C, it includes the product ADP and GlcNAc-1-P, which shows closed state on basis of MD simulations, implying that the protein maintains closed during the chemical reaction in catalysis. As to Model D, the protein does not contain GlcNAc, and the pocket becomes smaller after MD simulations. Moreover, the comparisons of Models A & B and Models A & D have been also performed, and we found that the lid motif has a closing tendency for Model B, and maintains the closed state for Model A and Model D. This means that the protein prefers to keep a closed state as the ligand stays at the active site. However, in experiment, Sato et al.10 obtained the crystal of complex of NahK•ATP in open state. Presumably, the GlcNAc is just ready to enter the protein or just release when the lid motif had no enough time to shut down its pocket during crystallization of the complex. In order to gain insight into such structural features of the active domain, the transportation of GlcNAc and related conformational changes of protein have been investigated and discussed (vide infra).

Figure 3. The overlap of initial structures and MD simulations structures for Models A, B, C, and D, as well as the comparison of Models A & B and Model A & D. The ellipse is “lid motif”, which refers to loop 26 ~ 32 and loop 314 ~ 321.

3.2. Identification of channels The conformations of Model A and Model B are extremely similar after 20 ns MD simulations, and the selected protein ID of Model A is 4OCK, which not only has better resolution and more complete amino acid sequence, but is also more close to the favorable conformation of NahK for its closed structure when GlcNAc and

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 18

Table 1. The location, number, and possibility of trajectories for GlcNAc delivery pathways on basis of RAMD MD simulations. The possibility was estimated by the ratio of Nx/Ntotal, where Nx is the number of trajectories of channel Px (x = a ~ h), and Ntotal is the number of total trajectories. Pathway

Location

Number

Possibility

Pa

Between residues of loop 4 and loop 6

25

39.0%

Pb

Between residues of loop 2 and loop 5

19

29.7%

Pc

Among residues of loop 2, loop 5, and loop 7

7

10.9%

Pd

Among residues of loop 1, loop 2, and loop 7

6

9.4%

Pe

Among Residues of loop 6

3

4.7%

Pf

Between residues of loop 6 and helix 2

2

3.1%

Pg

Between residues of helix 1 and helix 2

1

1.6%

Ph

Between residues of loop 2 and loop 3

1

1.6%

ATP stay at the active site. We thus analyze the equilibrium structure of the enzyme complex with GlcNAc and ATP on basis of Model A. As we can see in Figure S2, when the system arrives at a stable state, the GlcNAc breaks away from ATP. Here the GlcNAc is stabilized by Asp208, Asp243, Glu254, Arg306, and Tyr317 through eight hydrogen bonds. The ATP is localized by the coordinate bonds from two magnesium ions and hydrogen bond interactions from Asn33, Gln48, Tyr102, Lys210, and Asn212. It is worth noting that the Lys210 connects both GlcNAc and ATP directly and indirectly, which provides interactions with two molecules at the active site, and controls them do not to be far apart from each other. The conformation of the complex for NahK•ATP•GlcNAc will be used for the identification of transport channels in calculations for GlcNAc.

On basis of RAMD MD simulations, 64 trajectories are obtained and classified into eight possible pathways, which named as Pa, Pb, Pc, Pd, Pe, Pf, Pg, and Ph, as shown in Table 1 and Figure 4. Herein, Pa is located between residues of loop 4 and loop 6; Pb lies between residues of loop 2 and loop 5; Pc is seated among residues of loop 2, loop 5, and loop 7; Pd is situated among residues of loop 1, loop 2, and loop 7; Pe is located at residues of loop 6; Pf involves residues of loop 6 and helix 2; Pg lies between residues of helix 1 and helix 2; Ph is seated between residues of loop 2 and loop 3. The possibilities of these channels are 39.0%, 29.7%, 10.9%, 9.4%, 4.7%, 3.1%, 1.6% and 1.6% for Pa, Pb, Pc, Pd, Pe, Pf, Pg, and Ph, respectively. The most predominant channels, Pa, will be analyzed in detail in the following discussion. 3.3. The thermodynamic and dynamics character of GlcNAc transportation The transportation coordinate is defined as distance between CA of Arg196 and C5 of GlcNAc, which named as TC1 (see Figure S3), orientated to the door of Pa. Interestingly, the substrate is more likely transported from the entrance comprised with four loops of residues 26 ~ 32,

Figure 4. The key residual segments of possible channels for GlcNAc transportation on basis of RAMD MD simulations. (Loop 1 (red): 26-32; Loop 2 (orange): 51-57; Loop 3 (pink): 89-93; Loop 4 (blue): 109-115; Loop 5 (violet): 145-147; Loop 6 (green): 248-261; Loop 7 (yellow): 314-321; Helix 1 (nattierblue): 295-300; Helix 2 (celeste): 331-345)

Figure 5. The free energy profile of GlcNAc transportation along the TC1 .

ACS Paragon Plus Environment

Page 7 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Figure 6. The structures of protein and conformational changes of key residues along the TC1 in the GlcNAc transportation.

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

109 ~ 115, 248 ~ 261, and 314 ~ 321. The process of GlcNAc binding can be approximately classified into four stages on basis of the free energy changes and binding mode evolutions, which are shown in Figures 5 and 6. Moreover, in order to verify the convergence of PMF curves47, different time period and district are considered (see supporting information paragraph 1.1 and Figure S4). In the first stage (TC1 > 36.5 Å), the substrate hovers around the protein entrance, and it exposes the water environment without any binding force from the residues. Herein, the free energy has a weak fluctuation. In the second stage (34.2 Å < TC1 ≤ 36.5 Å), the hydrogen bond interaction has been observed, which drives the substrate access into the protein. For example, at 36.5 Å, the side chain of Asp255 provides one hydrogen bond with GlcNAc, which is benefit for the substrate recognizing the entrance of channel. Since the interaction is detected, the free energy is going down slightly at the point. However, the binding force is too small to fix the substrate stay at one location, and the other residues of Asp29, Asn212, and Glu254 as well as ATP molecule also form hydrogen bonds with GlcNAc at different time points (see Table S2). They can be also responsible for the door recognition of protein for the substrate, and guides free energy profile to show a decrease tendency with some fluctuation. In the third stage (29.7 Å < TC1 ≤ 34.2 Å), more interactions are detected than the second stage, which facilitate the substrate moving along the channel, and thus the free energy is apparently lower than the first stage. Herein, there are four flexible loops around the substrate, and thus the interactions are movable. The stage experiences three different periods (see supporting information 2). One feature should be noted that the GlcNAc experiences several distinct flips during the stage, which probably caused by the unstablility of the loops around it. Herein two to four hydrogen bonds can be observed, which come from different residues of Asp29, Asn112, Lys210, Arg246, Glu254, and Arg329, (see Figure 6 and Table S2), driving the free energy transient oscillation or steady decent. In the fourth stage (25.0 Å ≤ TC1 ≤ 29.7 Å), the substrate accesses into the active site step by step, and becomes much more stable by the assistance of more stable hydrogen bond interactions. With the binding force enhancement, the free energy is apparently going down to arrive the minimal point. Herein, Asp208, Lys210, Glu254, and Arg306 residues are involved in formation of hydrogen bonds with polarized groups of GlcNAc. As we can see in Figure S5, the binding modes are similar at different points, and more hydrogen bond interactions have been detected than the first three stages. The strong interactions are benefit for the stabilization of protein conformation and the substrate orientation.

Page 8 of 18

by around 0.8 and 2.5 kcal/mol, respectively. The fourth stage is of remarkable exothermicity, and the free energy descends about 6.5 kcal/mol. The tendency of energy changes is highly dependent with the quantity and quality of hydrogen bonds between residues and GlcNAc. The stronger the hydrogen bond and the more the number, the bigger the driving force of substrate binding and the faster the free energy falls. 3.4. The conformational changes of lid motif In experimental study, Sato et al.10 found that the protein experiences conformational changes from “open” to “closed” upon the ligand binding. In our present study, the thermodynamic and dynamics characters of the main channel for the GlcNAc transportation have been identified. The GlcNAc release process is primarily concerned with Loops 1, 4, 6, and 7 (see Figure S6), and thus the closest distances for the main chain of Loops 1 ~ 4, Loops 4 ~ 6, Loops 6 ~ 7, and Loops 7 ~ 8 along the transportation coordinate are summarized in Figure 7. It can be seen that the biggest distances of Loops 1 ~ 4, Loops 4 ~ 6, and Loops 1 ~ 7 appear in the third stage; while for the distance of Loops 6 ~ 7 maintains around one value, implying that the protein probably open its pocket at the third stage, which highly relevant to the distance of Loops 1 ~ 4, Loops 4 ~ 6, and Loops 1 ~ 7. Therefore, four representative snapshots along the reaction coordinate that covers four phases of the substrate binding are extracted to discuss the conformational changes of “lid motif” deeply, which also lend support to the prediction above (see supporting information 3 and Table S3).

Figure 7. The distances for lid motif and overlap of NahK for different stages along the TC1 in GlcNAc delivery (First stage: silver; Second stage: red; Third stage: cyan; Fourth stage: green).

The whole process of substrate binding is an obvious exothermic process, and the binding free energy is about 9.6 kcal/mol. It will provide a driving force for the following chemical reactions and product release. During the second and third stages, the free energy decreases slowly

ACS Paragon Plus Environment

Page 9 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Figure 8. The reaction scheme, free energy profile of GlcNAc-1-P production along RC2 as well as the structures of reactant, transition state, and product for Model E on basis of QM/MM MD simulations. In the reaction scheme, the residues involved in 2+ Pi and proton transfer are listed. The main coordination between Mg and ligands and the transferred groups are highlighted in dashed and blue.

Experimentally, Sato et al.10 put their “open” NahK•ATP complex crystal structure overlapped with “closed” NahK•AMPPNP•GlcNAc complex obtained by Wang et al.9 with N-domain and C-core domain align, and detected that the C-helical domain has an open tendency for NahK•ATP complex compared with NahK•AMPPNP•GlcNAc complex. Herein, we also perform an “overlap study” for the four key states above. As we can see in Figure 7, the conformational changes of Chelical domain are mainly related with loops 1 and 4. When the N-domain and C-core domain are overlapped, the conformation of the C-helical domain almost aligned for the first and second stages. Afterwards, a clockwise rotation has been observed for the C-helical domain in the third stage, and an open state has been detected for the lid motif. Then the C-helical domain recovered gradually in the fourth stage, which overlapped with the structure of first stage, showing a closed state. Therefore, the predicted conformational changes of lid motif between “open” and “closed” during the GlcNAc transportation agree with the experimental study, and the open-state structure of the NahK complex with ATP obtained by Sato et al.10 is probably located at the third stage. 3.5. Chemical reaction mechanism for GlcNAc-1-P production The generation of GlcNAc-1-P is probably related with the transfer of phosphate group and proton, and thus the trend of GlcNAc-1-P production is shortening of H1-O3 and O2-Pγ distance and stretching of H1-O2 and O4-Pγ. Herein two reaction coordinates are tested, which defined as dH1-O2 + dO4-Pγ- dH1-O3 - dO2-Pγ (RC1) and – dH1-O3 – dO2-Pγ (RC2). It can be seen from Figure S7, both of them show

almost the same relative energy profiles. Considering the computational cost, the RC2 has been selected as reaction coordinate in free energy calculations. The trajectories of different time periods are tested to guarantee the PMF convergence (see supporting information 1.2 and Figure S8). Herein, the profiles of different time periods of 16 ~ 19, 17 ~ 20, 16 ~ 20 ps overlap completely, which means that the convergence of QM/MM MD simulations for the PMF curves is reliable. The potential of mean forces (PMFs) and all of the corresponding important states are depicted in Figure 8. The PMF profiles show that the GlcNAc deprotonation happens synchronously with the phosphate group transfer from ATP to GlcNAc, since only one transition state has been discovered. Due to its low free energy barrier (6.9 kcal/mol) as well as exothermicity, the reaction is kinetically and thermodynamically reasonable. Moreover, the barrier is lower than the value of 16.2 kcal/mol estimated by the transition state theory from experimental data (kcat is 9.9 ± 0.3 s-1, which measured at 4.0 mM MgCl2), indicating that the rate-limiting step probably is arisen from the other stages. As we can see in Figure 8, at the reactant state, strong ion-pair bonding interactions are detected between the NH3+ group of the protonated Lys210 and phosphate group of ATP, and a hydrogen bond is found between Asn213 and phosphate group; moreover, another hydrogen bond can be identified between Asp208 and GlcNAc. The interactions are helpful in localizing two substrates, which are favorable for the reaction occurring. From reactant to transition state, the distance of O2 and Pγ is shortening from 2.99 ± 0.08 Å to 2.17 ± 0.08 Å, together with the elongating of O4 and Pγ from 1.84 ± 0.07 Å to 2.55 ± 0.16 Å. However, at the same time, the changes of the distance between O2 and H1 are small, which just

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

increase from 1.00 ± 0.02 Å to 1.04 ± 0.03 Å. The phenomenon implies that the phosphate group transfers in preference to proton. At the product state, the bond of O2- Pγ and O3-H1 are formed entirely, whose bond lengths are 1.70 ± 0.03 Å and 1.01 ± 0.02 Å, respectively. In order to verify the transfer sequence of phosphate group and proton, the entire snapshots along the reaction coordinate are gathered to analyze. As shown in Figure S9, the crossing point of O3-H1 and O2-H1 localizes at about - 3.15 Å, while those of O2-Pγ and O4-Pγ stay at ~ - 3.70 Å of RC2. In addition, both O3-H1 and O2-Pγ are bonding at the ~ 2.75 Å of RC2. Accordingly, the transfer of phosphate group takes precedence over GlcNAc deprotonation, while both transfer processes are completed almost at the same time. It should be noted that the O2-H1 distance distribution has two distinct values around the 3.0 Å of RC2, and thus the snapshots are collected from -3.01 Å to 2.99 Å along the reaction coordinate (see Figure S10). The vast majority of O3-H1 distance locates at 1.0 ~ 1.1 Å, which obviously shorter than that of O2-H1, showing that the O3-H1 bond is formed or has a bonding trend. Therefore, if the orientation of O3-H1 bond changes, the distance of O2-H1 will be greatly influenced, and thus large fluctuation will be observed for distance of O3-H1. Moreover, another interesting thing should be noted. Sato et al.10 predicted that the possible role of Asp208 as a positioning residue for correcting the location of anomeric hydroxyl group of substrate rather than facilitating the deprotonation of substrate. Here MM MD and QM/MM MD simulations have been performed on the Asp208Ala mutated system (see supporting information 4). Based on the stabilized structure after 20 ns molecular dynamics simulations, we found that the direction of anomeric hydroxyl group is changed, which towards to the Pi group of ATP (see Figure S11). Clearly, the Asp208 can stabilize the position of substrate for nucleophilic attack. As to its catalytic role, the QM/MM study shows that the energy barrier for the Asp208Ala mutated system is about 50 kcal/mol, which is much higher than that of wild system by 41 kcal/mol (see Figure S12). It means that the Asp208Ala mutation may lower the catalytic activity of the enzyme to some extent; while the reaction might be also happened, and ATP supports oxygen for nucleophilic attack. As to the coordination of the two magnesium ions, six coordination bonds are detected during the whole reaction. As we can see in Figure S13, the first Mg2+ bridged the α-γ phosphates of ATP coordinates with OPα (the terminal oxygen of Pα), OPγ (the terminal oxygen of Pγ), O4 (the bridging oxygen of Pβ-Pγ), Asn213, Asp228, and one water molecule. Herein, we mainly focus on the dative bonds of Mg1-OPα, Mg1-OPγ, and Mg1-O4, which are involved in the chemical reaction. The coordination bond of Mg1-O4 is apparently weaker than the others at the reactant state. As reaction happens, the bridging oxygen of Pβ-Pγ evolves into terminal oxygen of Pβ gradually, and the distance of the coordination bond decreases, which means that the coordination is strengthened. Meanwhile, the coordination bond of Mg1-OPα is gradually weakened,

Page 10 of 18

which is probably caused by the coordination of Mg1-O4. With regard to the dative bond of Mg1-OPγ, as the phos phate group transfer, its distance is stretching from 2.04 ± 0.06 Å to 2.20 ± 0.11 Å, which implies that the coordination is weakening. For the second Mg2+ bridged the β-γ phosphates of ATP, it coordinates with OPβ and OPγ of ATP, two oxygen atoms of carboxyl group in Asp228, and two waters. Herein, the Mg2-OPβ and Mg2-OPγ are relative with the phosphate group transfer, while both of them vary slightly. The distance of Mg2-OPβ is shortened from 2.09 ± 0.08 Å to 2.04 ± 0.06 Å, which means that the coordination becomes weak feebly. For Mg2-OPγ, the coordination interaction is enhanced from the reactant state to the transition state, while it is then weakened and its distance at the reactant and product states is almost the same. Here we can see that the Mg1 is more sensitive to the phosphate group transfer than Mg2. Moreover, our calculated NahK•ADP•GlcNAc-1-P and NahK•ATP•GlcNAc complex with two magnesium ions are used for comparing with the crystal structure of NahK•ADP•GlcNAc-1-P complex with one magnesium ion and NahK•ATP complex with two magnesium ions obtained by experiment9-10 (see Table S4 and Figure S14), we find that the position of the Mg2+ and the coordination in our built model are perfectly consistent with those of crystal structure containing the Mg2+, which means that our QM/MM model is reasonable, and the deficient Mg2+ probably refers to the Mg2. Sato et al.10 found that the secondary magnesium ion is important for catalysis of NahK, and thus we will analyze the role of it in detail in the next discussion. 3.6. The role of the secondary Mg2+ We mainly discuss the role of the secondary Mg2+ that bridging the β-γ phosphates on basis of Model F. For theproduction of GlcNAc-1-P, the reaction coordinate is defined as – dH1-O3 – dO2-Pγ (RC2’), which is the same as RC2. The trajectories during different time periods are collected to analyze the convergence of PMF (see supporting information 1.3 and Figure S15). The free energy profile and corresponding key states of reactant, transition state, intermediate, and product are listed in Figure 9. We can see that the catalytic reaction can be classified into two steps: (i) The phosphate group transfers from ATP to GlcNAc, together with the deprotonation of GlcNAc mediated with Asp208 and Asp228; (ii) the proton transfers from Asp228 to Asp208. At the reactant state, the Mg2+ is tetrahedrally coordinated with the Pα, Pγ, Asn213, and Asp228, which is different from Model E that is hexa-coordinate. Here the distance of Pγ and O2 is 3.87 ± 0.09 Å, and a hydrogen bond is detected between O2 of GlcNAc and carboxyl of Asp208 with the distance of 1.64 ± 0.09 Å. At TS1c’ state, the O2 ••• Pγ and H1 ••• O3 are 2.35 ± 0.16 Å and 1.68 ± 0.16 Å, respectively. It likes that the phosphate group transfers prior than the proton delivery. Then the phosphate group is still moving toward to GlcNAc, and the proton departs from GlcNAc, which shifts to Asp228 mediated by Asp208.

ACS Paragon Plus Environment

Page 11 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Figure 9. The reaction scheme, free energy profile of GlcNAc-1-P production along RC2’ as well as the structures of reactant, intermediate, transition state, and product for Model F on basis of QM/MM MD simulations. In the reaction scheme, the resi2+ dues involved in Pi and proton transfer are listed. The main coordination between Mg and ligands and the transferred groups are highlighted in dashed and blue.

At IM1c’ state, the GlcNAc-1-P is generated, and the distance of O2 ••• Pγ stays at 1.68 ± 0.03 Å. Herein the Mg2+ is penta-coordinated with the terminal oxygen of GlcNAc-1P, OPα and OPγ of ADP, Asn213, and Asp228 (see Figure S16). Afterwards, from IM1c’ to Pc’ state, the proton transfers from Asp228 to Asp208, and the distance of proton and carboxyl of Asp228 decreases from 1.62 ± 0.03 Å to 1.01 ± 0.03 Å. The coordination is less changed, except for the reinforce of Mg-OAsp228, which probably caused by the proton transfer. In order to explore this step in depth, the corresponding trajectories are employed to further analysis. The last 5 ps snapshots of each window are selected, which are shown in Figure S17. From - 5.5 Å to – 3.7 Å along the reaction coordinate, the distance between O2 of GlcNAc and Pγ of ATP decreases, while the H1 hovers between O2 and O3, and not form interaction with Asp228. It means that the phosphate group transfers before GlcNAc deprotonation at the initial stage. From – 3.7 Å to – 3.3 Å, the distance of O2 and Pγ fluctuates between 1.6 and 2.0 Å at first, and then stabilizes at 1.6 Å gradually. The distance of O2 ••• H1 stretches from 1.0 Å to 3.5 Å hastily, together with shortening and bonding of OAsp228 ••• H1. It implies that phosphate and proton transfer almost completely at the same time. From -3.3 Å to – 2.7 Å, the distance of O2-Pγ almost not change, while the pro-

ton transfer from Asp228 to Asp208 gradually arrives the stable state. Herein, two distinct values for O2-H1 and OAsp228-H1 occur between – 3.7 Å to – 3.5 Å along the reaction coordinate, and thus the distance of O2-H1, O3-H1, and OAsp228-H1 as well as the charges of O2, O3, and OAsp228 during the region are collected for analysis (see supporting information 5 and Figure S18). We find that the phenomenon probably comes from the random orientation of Asp208, which plays a mediated role in the proton transfer. Herein, one question has been arisen: why does not the proton transfer to Asp208 directly and instead wanders between Asp208 and Asp228? To solve it, the statistic charges of O3 and OAsp228 of the related snapshots are collected in Figure S19, along the RC2’ from – 3.7 Å to – 3.3 Å. The statistic charges of O3 mainly locate at around -0.55 a.u., which are more negative than those of OAsp228 by 0.1 a.u., and the charges of majority snapshots from – 0.7 a.u. to – 0.1 a.u. overlaps, which imply that the O3 has stronger attraction toward the proton, while the proton affinity of Asp228 should not be neglected. To obtain more detailed information, the charges of relevant snapshots along the reaction co-ordinate are also analyzed. The overlap region mainly locates from - 3.7 Å to - 3.3 Å, sug-

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

gesting the charges of both OAsp228 and O3 are comparative and the proton prefers to wander between O3 and OAsp228. This also creates the condition for the proton transfer to Asp228. However, in the second step, the protonated Asp228 easily gives the proton to Asp208. From 3.2 Å to – 2.8 Å, the charges of O3 are more negative than those of OAsp228, and the overlap region is obviously reduced, which also accounts for the phenomenon. With regard to the free energy changes, we can see that if we only consider one magnesium ion, the chemical reaction is very different with that of two magnesium ions. For the catalytic mechanism, the reaction follows a stepwise mechanism, while for Model E, it shows a concerted

Page 12 of 18

mechanism. As to the free energy barrier, the ratelimiting step is about 16.5 kcal/mol, which is higher than that of Model E by 9.6 kcal/mol. Therefore, the 200 snapshots at - 4.7 Å (Model E) and - 5.5 Å (Mode F) of the reaction coordinate have been considered, and the relative charges are collected. As we can see in Figure 10, the charges of O4 of ATP are more negative than O2 of GlcNAc by 0.60 a.u. for Model E, while for Model F, the charges of O4 are positive than O2 by 0.2 a.u. It implies that the transfer will easily occur for model E compared with model F, which agrees with the free energy profiles. Moreover, the whole phosphate group is involved in the transfer process, and thus the charges come from the phosphate group will be considered. The charges of phosphate group for Model E maintain at - 1.14 ± 0.13 a.u., while for Model F, the charges are -1.42 ± 0.09 a.u.. Besides that, the charges of O2 for Model E are more positive than Mode F, and the charges of O4 are more negative than those of Model F. The more negative charges for oxygen atom are, the stronger repulsive interaction is. Therefore, the charge distribution is favorable for the phosphate group transfer to some extent towards for Model E. Moreover, in the reactant state, the distance of Pγ and O2 is 2.99 ± 0.08 Å, shorter than that of Model F by 0.88 Å. The difference of bond distance probably comes from the charge distribution, which provides a better condition for the phosphate group transport. It thus can be concluded that the secondary magnesium ion that bridging the β-γ phosphates is necessary in the chemical catalytic process. 3.7. The protonation of GlcNAc-1-P It is worth noting that the distance between H1 of protonated Asp208 and Pγ of GlcNAc-1-P is 1.68 ± 0.12 Å. If the proton transfers from Asp208 to GlcNAc-1-P, the binding of one magnesium ion may be relaxed. Hence, in order to obtain the thermodynamic and dynamic characterizations, the distance stretching of OPγ-Mg1 and H1-O3 has been defined as the reaction coordinate (RC3= dOPγ-Mg1 + dH1-O3). Based on the QM/MM MD simulations combined with the umbrella sampling technique, the potential of mean forces is tested during different time periods to ensure the convergence (see supporting information 1.4 and Figure S20). The last PMF curves as well as the corresponding key states are obtained, which are listed in Figure 11.

Figure 10. The statistic charges of main atoms and groups for reactant state in model E and F.

As Figure 11 shows, the discrepancy among the free energy profiles from different time periods are very small, implying that reliable convergence of PMF curves are achieved. Clearly, the free energy barrier for this step is very small, and the reaction is slightly endothermic, suggesting that the proton might easily transfer between Asp208 and GlcNAc-1-P. Moreover, we note that the fluctuation of curves is large along the reaction coordinate, which implies that the structures of the snapshots during the proton transfer are probably unstable.

ACS Paragon Plus Environment

Page 13 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

Figure 11. The reaction scheme, the free energy profile of GlcNAc-1-P protonation along RC3 as well as the structures of key states on basis of QM/MM MD simulations. In the reaction scheme, the residues involved in Pi and proton transfer are listed. 2+ The main coordination between Mg and ligands and the transferred groups are highlighted in dashed and blue.

At the GlcNAc-1-P deprotonated state, the distances of OPγ-H1 and OPγ-Mg1 are 1.68 ± 0.12 and 2.17 ± 0.03 Å, respectively. As the reaction proceeds, the H1 transfers from the protonated Asp208 to GlcNAc-1-P, and the OPγ gradually gets rid of the magnesium ion. At the product state, the distances of OPγ-H1 and OPγ-Mg1 are 1.06 ± 0.04 and 2.48 ± 0.10 Å, respectively. For the sake of identifying the stability of various key states and the specific changes of bond length along the reaction coordinate, the snapshots for all windows are collected to make a discussion. As shown in Figure S21, the proton transfers from Asp208 to OPγ synergistically happened with OPγ breaking away from the magnesium ion. The distance of H1 and OPγ de creases gradually, and its fluctuation is also weakened. At the same time, the distance of O3 and H1 increases, and the fluctuation of bond length is reinforced. It shows that the proton wanders between OPγ and O3, implying that the protonation of GlcNAc-1-P may be not stable, which agrees with the endothermic phenomenon. As to the distance of OPγ and magnesium ion, it stretches piece by piece, while its fluctuation enhances gradually along the reaction coordinate, especially after RC value of 3.7 Å. This suggests that with the OPγ breaking away from magnesium ion, the corresponding coordination bond becomes instability, which also indicates that the step is probably endothermic. We can see that the free energy profile reaches convergence, while the undulations are continual, which may be also mainly caused by the distance fluctuations of OPγ-H1, O3-H1, and OPγ-Mg1. With regard to the two magnesium ions, six dative bonds have been observed during the proton transfer process, which are listed in Figure S22. For the Mg2, the length of coordination bond is almost unchanged along the reaction coordinate; while for the Mg1, the distance of OPγ-Mg1 gradually increases as the proton transfer from Asp208 to OPγ.

As can be seen from the above discussion, the protonation of GlcNAc-1-P is benefit for it flees away from the binding of magnesium ion, while the reaction is not stable, and a weak interaction is also existed between protonated oxygen and Mg2+. Probably, the effect of protonation of GlcNAc-1-P is not obvious enough for GlcNAc-1-P broken away from Mg2+ binding. Interestingly, several water molecules have been detected around the magnesium ions, constructing a favorable environment for the coordination with Mg2+ instead of GlcNAc-1-P (see Figure S23). 3.8 Water-assistant protonated GlcNAc-1-P release Herein, on basis of Phc state in QM/MM MD simulations, two water molecules, about 4 ~ 5 Å away from the two magnesium ions have been concerned , which are selected as the candidates of GlcNAc-1-P to coordinate with Mg2+. The shortening of the distances Mg1-OWAT2 (RC4) and Mg2-OWAT5 (RC5) are defined as reaction coordinates, which listed in Figure 12. After a series of MM MD simulations (see supporting information 6 and Figure S24), the 2D-PMF curve is obtained (see Figure 13), and three possible mechanisms have been detected. (I) The Mg1 and OWAT2 form the dative bond prior to that of Mg2 and OWAT5, which overcomes two free energy barriers of 9.4 and 7.3 kcal/mol in binding sequence; (II) The coordination bond of Mg2 and OWAT5 are formed before Mg1 and OWAT2, and the free energy barriers are 8.3 and 12.9 kcal/mol, respectively; (III) The Mg1-OWAT2 and Mg2OWAT5 form coordination bonds simultaneously. The most favorable mechanism will be discussed in detail, and the coordination structures of two magnesium ions for reactant, intermediate, transition state, and product are shown in Figure 12. At Rw state (RC4 = 5.2 Å, RC5 = 4.2 Å), six coordination bonds are detected for both Mg2+. For Mg1, it coordinates with protonated OPγ of Glc-

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

NAc-1-P, OAsn213, OAsp228, WAT1, OPα, and OPβ. Two distances of Mg1-OWAT2 and Mg1-OPγ involved in the process of exchange for protonated GlcNAc-1-P and WAT2 are 5.20 ± 0.06 and 2.28 ± 0.12, respectively. As to Mg2, six coordination interactions from deprotonated OPγ of GlcNAc-1-P, OPβ of ADP, two oxygen atoms of Asp228, and two water molecules are found. We mainly focus on the two distances of Mg2-OWAT5 and Mg2-OPγ, which are 4.20 ± 0.03 and 1.91 ± 0.50 Å, respectively. As exchange progresses, the WAT2 closes to Mg1, and the protonated OPγ of GlcNAc-1-P flees away from the binding of Mg1 simultaneously. Here the coordination of Mg2 is not changed notably. At IM1w state (RC4 = 2.0 Å, RC5 = 4.2 Å), the distances of Mg1- OWAT2 and Mg1- OPγ are 2.01 ± 0.03 and 4.45 ± 0.14 Å, which means that the WAT2 completely replaces GlcNAc-1-P to coordinate with Mg1. For Mg2, it

Page 14 of 18

still maintains six coordination structure. However, one interesting thing happened. The oxygen atom of Asp208 coordinates with Mg2 instead of WAT4 with distance of 1.97 ± 0.07 Å. The phenomenon probably occurs by the rotation of GlcNAc-1-P leaving away from the Mg1 completely. After that, the WAT5 starts closing to Mg2. We found that the distance between deprotonated OPγ of GlcNAc-1-P and Mg2 maintains around 1.89 Å during the process, while distance of OAsp208-Mg2 is stretching as the binding of Mg1-OWAT5. At Pw state (RC4 = 2.0 Å, RC5 = 2.0 Å), the WAT5 instead of Asp208 to coordinate with Mg2 completely with the distance of 2.00 ± 0.03 Å, and the distance of Mg2-OPγ is 1.89 ± 0.04 Å. It can be seen that from IM1w to Pw state, it is an exchange process of WAT5 and Asp208, which means that the coordination bond of deprotonated oxygen atom and Mg2 is so strong that probably needs higher free energy to break it. For the second possible mechanism, the first free energy barrier mainly comes from the exchange of WAT5 and WAT4 that binds with Mg2, and the second barrier chiefly derived from the substitution of GlcNAc-1-P by WAT2 that coordinates with Mg1. The different free energy barriers for the mechanisms I and II are probably caused by the distinct coordination interactions of both Mg2+ ions in IMw1 and IMw2 (see Figure S25). For the third possible mechanism, the free energy barrier mainly comes from the unstable penta-coordination structure of both magnesium ions, which is listed in Figure S26.

Figure 13. The two-dimensional free energy profile for GlcNAc-1-P cleavages from the magnesium ion.

Figure 12. The length of coordination bonds for Mg1 and Mg2 in mechanism I. (The RC4 and RC5 are highlighted in 2+ red. The distances between GlcNAc-1-P and Mg are highlighted in blue)

Therefore, the coordination of Mg2 and deprotonated oxygen atom of GlcNAc-1-P requires higher free energy barrier to assist its cleavage. Moreover, as to the other product ADP, it has three coordination bonds with two Mg2+, which probably more difficult than GlcNAc-1-P for breaking away from the binding of Mg2+. It means that the Pi-Pi breaking in the chemical reaction probably not the rate-limiting step during the whole catalytic process. Moreover, we detected that there is a hydrogen bond between protonated oxygen atom of GlcNAc-1-P and carboxyl group of Asp208, and thus Asp208 possibly as “pro-

ACS Paragon Plus Environment

Page 15 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

ton-shuffle” to mediate the proton transfer to the deprotonated oxygen atom of GlcNAc-1-P binding with Mg2, which is benefit for the product flee away from Mg2. Besides that, there is also another hydrogen bond between His31 and another deprotonated oxygen atom of GlcNAc1-P, and thus the transfer of proton may alter the charge distribution of GlcNAc-1-P, which is probably useful for the cleavage of GlcNAc-1-P (see Figure S27). The completely cleavage of GlcNAc-1-P and ADP as well as the kinetic characteristics of their release entirely from

the protein will need more studies both experimentally and theoretically. For example, Markov State Model (MSM) is a powerful tool for accessing long time scale dynamics from many short MD simulations in substrate transportation48-49. In our system, the microstates can be generated by clustering the seeds from unbiased MD simulations on basis of the steered MD (SMD), targeted MD (TMD), and RAMD MD etc., which covers the whole release process, and then macrostates can be identified and capture long time scale dynamics to analyze the kinetics of GlcNAc-1-P or ADP flee away from the protein.

Figure 14. The whole free energy profile for GlcNAc transportation, GlcNAc-1-P production and protonation, and water assisted 2+ GlcNAc-1-P cleavage from Mg on basis of MM MD and QM/MM MD simulations.

4.

CONCLUSION

Based on the extensive MM MD and QM(DFT)/MM MD simulations, two magnesium ion dependent enzymatic catalysis of NahK and its substrate delivery have been studied. NahK is the first example of an anomeric kinase acting on gluco-type substrate, and its two magnesium ions are structurally similar with some protein kinases and aminoglycoside phosphotransferases. In the present study, an initial picture of the catalytic process has been built (see Figure 14), and several interesting facts have been revealed. Firstly, in chemical reaction of the GlcNAc-1-P production, a concerted mechanism has been identified for the delivery of phosphate group and proton, while the phosphate group transfers from ATP to GlcNAc is prior to the delivery of proton from GlcNAc to Asp208. This is an exothermic reaction with the free energy barrier of 7.0 kcal/mol, which is lower than that of GlcNAc-1-P cleavage from magnesium ions by 2.4 kcal/mol, suggesting that the cleavage of Pi-Pi bond is possibly not the rate-limiting step in the enzymatic process. This is quite different from the other active site structurally similar kinases, whose Pi-Pi bond activation almost determines the whole catalytic efficiency50-51. The two magnesium ions always maintain six coordination with GlcNAc (GlcNAc-1-P), ATP (ADP), and Asp228, Asn213, and water molecules, although small fluctuation has been observed. Moreover, the protonation mechanism of GlcNAc-1-P has been discovered to play a benefit role in getting rid of one magnesium ion binding, and combined with the assistant of water, the protonated oxygen atom in GlcNAc-1-P can break away from Mg1 completely. The second significant

thing happens during the substrate transport process. Eight possible pathways have been identified. With regard to the most favorable channel, the substrate binding is remarkable exothermic, which is mightily dependent on the quantity and quality of hydrogen bond network around the ligand. The stronger the hydrogen bonds, the bigger the driving force for the GlcNAc binding, and the corresponding exothermic trend is more obvious. Four stages have been classified, based on the binding mode and free energy changes. The channel is related with four flexible loops, and thus different residues may assist the entrance of substrate at different stages. Moreover, an open-closed conformational change has been identified upon the ligand binding, which occurs in the third stage. Another interesting finding should owe to the role of secondary Mg2+ that bridging the β-γ phosphates of ATP. The deficiency of the magnesium ion will cause different mechanism and higher free energy barrier, which reduces the enzymatic efficiency visibly, implying that the Mg2 play an active role in the catalysis. Moreover, the whole free energy profiles for different steps, including GlcNAc transportation, GlcNAc-1-P production and protonation as well as the water-assisted GlcNAc-1-P cleavage from Mg2+, show that the present highest barrier of 9.0 kcal/mol derived from the broken of GlcNAc-1-P and Mg2+ is still lower than the experimental estimation of 16.2 kcal/mol. Two possible reasons are raised as follows: the theoretical modeling may have its limitation to some extent, and the following processes, such as GlcNAc-1-P and ADP fleeing away from both two Mg2+ and delivery to the outside of protein play an extremely important role in the

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

enzymatic catalysis, which probably overcome higher free energy barrier. The tentative test of water-assistant GlcNAc-1-P cleavage from binding of magnesium ions also proves that the deprotonated oxygen atom of GlcNAc-1-P getting rid of magnesium ions is more difficult than that of protonated oxygen atom, which lends support to this speculation. The new findings will supply a basis for comprehending other members of anomeric kinase, which is also benefit for understanding two magnesium ions structurally similar protein kinases and aminoglycoside phosphotransferases. The present work is thus useful for the further study of GNB/LNB pathway in bifidobacteria, and provides significant clues for biosynthesis of oligosaccharides and sugar 1-phosphates.

ASSOCIATED CONTENT Supporting Information Convergence test of PMF curves for GlcNAc transportation, GlcNAc-1-P production, GlcNAc-1-P protonation; Three different periods during the third stage in GlcNAc transportation; Analysis of four representative snapshots covers four phases of GlcNAc transportation; Role of Asp208; Explains of two distinct values for O2-H1 and OAsp228-H1 in GlcNAc-1-P production on basis of Model E; Calculated methods for Glc2+ NAc-1-P cleavages from Mg ; The original pdb ID, ligands, number of magnesium ions, MD simulation time, and conformational changes of lid motif for four models; The number of hydrogen bonds formed between residues (ATP) and GlcNAc as well as percent for each window by using the last 10 ns trajectories in MD simulations; The distances of the closest atoms on the main chains for loop 1 ~ 4, loop 4 ~ 6, loop 6 ~ 7, and loop 1 ~ 7 obtained by using the last 10 ns trajectories in MD simulations; The crystal structures obtained by experiment; The root mean square deviation (RMSD) of the protein backbone for four models; The equilibrium structure of NahK•GlcNAc•ATP and its active site on basis of Model A; The defined reaction coordinate (TC1) for GlcNAc transportation; The structures of protein and conformational changes of key residues along the TC1 in the fourth stage of GlcNAc transportation; The location of four main loops that play an important role in GlcNAc transportation and the surface of the protein; Relative energy along two tested reaction coordinates in GlcNAc production; The statistics of main bonds that involved in Glc-NAc-1-P production and protonation along RC2 for Model E (F); The QM/MM model for Asp208Ala system in GlcNAc-1-P production; The coordination of two magnesium ions for reactant, transition state, and product in GlcNAc-1-P production and protonation for Model E (F) on basis of QM/MM MD simulations. The comparison for QM/MM models and experimental crystals in coordination of Mg1 and Mg2; The statistic charges of O3 and OAsp228 along the RC2’ from – 3.7 Å to – 3.3 Å as well as the statistic charges of O3 and OAsp228 along the RC2’; The twodimensional free energy profile for GlcNAc-1-P cleavages from magnesium ion during different time periods; The length of coordinated bonds for Mg1 and Mg2 in mechanism II and III; The hydrogen bonds formed between GlcNAc-1-P and residues of Asp208 and His31 at Pw state.

Page 16 of 18

AUTHOR INFORMATION Corresponding Author * E-mail: [email protected] (C.W) * E-mail: [email protected] (Z.C)

Notes The authors declare no competing financial interest.

Present Addresses †

The Key Laboratory of Natural Medicine and ImmunoEngineering, Henan University, Kaifeng 475004, People’s Republic of China ‡ State Key Laboratory of Physical Chemistry of Solid Surfaces and Fujian Provincial Key Laboratory of Theoretical and Computational Chemistry, College of Chemistry and Chemical Engineering, Xiamen University, Xiamen 360015, People’s Republic of China § College of Chemistry and Chemical Engineering, Henan University, Kaifeng 475004, People’s Republic of China

Author Contributions All authors have given approval to the final version of the manuscript.

ACKNOWLEDGMENT This work is supported by the National Science Foundation of China (21603057, 21373164, 21673185,), Project funded by China Postdoctoral Science Foundation, Project funded by Henan Postdoctoral Science Foundation, and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under Grant No.U1501501. We also thank Prof. Ruibo Wu at Sun Yat-sen University and Dr. Shenglong Wang at New York University for their assistance in computation. We also thank the National Supercomputing Center in Changsha (NSCCChangsha) and Guangzhou (NSCC-GZ) for providing the computational resources. This work was partially supported by State Key Laboratory of Physical Chemistry of Solid Surfaces (Xiamen University) (201618).

REFERENCES 1. Picard, C.; Fioramonti, J.; Francois, A.; Robinson, T.; Neant, F.; Matuchansky, C. Bifidobacteria as Probiotic Agents Physiological Effects and Clinical Benefits. Aliment. Pharmacol. Ther. 2005, 22, 495-512. 2. Kitaoka, M. Bifidobacterial Enzymes Involved in the Metabolism of Human Milk Oligosaccharides. Adv. Nutr. 2012, 3, 422S. 3. Kitaoka, M.; Tian, J.; Nishimoto, M. Novel Putative Galactose Operon Involving Lacto-N-Biose Phosphorylase in Bifidobacterium Longum. Appl. Environ. Microbiol. 2005, 71, 3158-3162. 4. Fushinobu, S. Unique Sugar Metabolic Pathways of Bifidobacteria. Biosci., Biotechnol., Biochem. 2010, 74, 2374-2384. 5. Nishimoto, M.; Kitaoka, M. Identification of NAcetylhexosamine 1-Kinase in the Complete Lacto-N-Biose

ACS Paragon Plus Environment

Page 17 of 18 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Catalysis

I/Galacto-N-Biose Metabolic Pathway in Bifidobacterium Longum. Appl. Environ. Microbiol. 2007, 73, 6444-6449. 6. Li, Y.; Yu, H.; Chen, Y.; Lau, K.; Cai, L.; Cao, H.; Tiwari, V. K.; Qu, J.; Thon, V.; Wang, P. G.; Chen, X. Substrate Promiscuity of N-Acetylhexosamine 1-Kinases. Molecules 2011, 16, 6396-6407. 7. Liu, Y.; Nishimoto, M.; Kitaoka, M. Facile Enzymatic Synthesis of Sugar 1-Phosphates as Substrates for Phosphorylases Using Anomeric Kinases. Carbohydr. Res. 2015, 401, 1-4. 8. Cai, L.; Guan, W.; Kitaoka, M.; Shen, J.; Xia, C.; Chen, W.; Wang, P. G. A Chemoenzymatic Route to N-Acetylglucosamine1-Phosphate Analogues: Substrate Specificity Investigations of NAcetylhexosamine 1-Kinase. Chem. Commun. 2009, 40, 2944. 9. Wang, K. C.; Lyu, S. Y.; Liu, Y. C.; Chang, C. Y.; Wu, C. J.; Li, T. L. Insights into the Binding Specificity and Catalytic Mechanism of N-Acetylhexosamine 1-Phosphate Kinases through Multiple Reaction Complexes. Acta Crystallogr. Sect. D. Biol. Crystallogr. 2014, 70, 1401-1410. 10. Sato, M.; Arakawa, T.; Nam, Y. W.; Nishimoto, M.; Kitaoka, M.; Fushinobu, S. Open–Close Structural Change upon Ligand Binding and Two Magnesium Ions Required for the Catalysis of N -Acetylhexosamine 1-Kinase. Biochim. Biophys. Acta, Proteins Proteomics 2015, 1854, 333-340. 11. Adams, J. A. Kinetic and Catalytic Mechanisms of Protein Kinases. Chem. Rev. 2001, 101, 2271. 12. Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: Automating pK Prediction and the Preparation of Biomolecular Structures for Atomistic Molecular Modeling and Simulations. Nucleic Acids Res. 2012, 40, W537-W541. 13. Myers, J.; Grothaus, G.; Narayanan, S.; Onufriev, A. A Simple Clustering Algorithm Can Be Accurate Enough for Use in Calculations of pKs in Macromolecules. Proteins: Struct., Funct., Bioinf. 2006, 63, 928-938. 14. Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A. H++: A Server for Estimating pKas and Adding Missing Hydrogens to Macromolecules. Nucleic Acids Res. 2005, 33, W368-W371. 15. Bashford, D.; Karplus, M. Pka's of Ionizable Groups in Proteins: Atomic Detail from a Continuum Electrostatic Model. Biochemistry 1990, 29, 10219-10225. 16. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935. 17. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157-1174. 18. Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The Resp Model. J. Phys. Chem. 1993, 97, 10269-10280. 19. Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. Gaussian 09; Gaussian Inc., Wallingford, CT, 2010. 20. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. 21. Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Bioinf. 2006, 65, 712-725. 22. Li, P.; Roberts, B. P.; Chakravorty, D. K.; Jr, K. M. M. Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent. J. Chem. Theory Comput. 2013, 9, 2733-2748.

23. Case, D.; Darden, T.; Cheatham III, T. E.; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Walker, R.; Zhang, W.; Merz, K. AMBHER 12; University of California, San Francisco, CA, 2012. 24. Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N⋅ Log (N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. 25. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593. 26. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327-341. 27. Lüdemann, S. K.; Lounnas, V.; Wade, R. C. How Do Substrates Enter and Products Exit the Buried Active Site of Cytochrome P450cam? 1. Random Expulsion Molecular Dynamics Investigation of Ligand Access Channels and Mechanisms. J. Mol. Biol. 2000, 303, 797-811. 28. Vashisth, H.; Abrams, C. F. Ligand Escape Pathways and (Un) Binding Free Energy Calculations for the Hexameric InsulinPhenol Complex. Biophys. J. 2008, 95, 4193-4204. 29. Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187-199. 30. Souaille, M.; Roux, B. t. Extension to the Weighted Histogram Analysis Method: Combining Umbrella Sampling with Free Energy Calculations. Comput. Phys. Commun. 2001, 135, 40-57. 31. Rosenbergl, J. M. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011-1021. 32. Becke, A. D. A New Mixing of Hartree–Fock and Local Density-Functional Theories. J. Chem. Phys. 1993, 98, 1372-1377. 33. Lee, C.; Yang, W.; Parr, R. G. Lee, C., Yang, W., Parr, R. G., Development of the Colle-Salvetti Correlation Energy Formula into a Functional of the Electron Density, Phys. Rev. B 37, 785 (1988). Phys. Rev. B: Condens. Matter 1988, 37, 785-789. 34. Zhao, Y.; Chen, N.; Wu, R.; Cao, Z. A QM/MM MD Study of the Ph-Dependent Ring-Opening Catalysis and Lid Motif Flexibility in Glucosamine 6-Phosphate Deaminase. Phys. Chem. Chem. Phys. 2014, 16, 18406-18417. 35. Zhao, Y.; Chen, N.; Mo, Y.; Cao, Z. A Full Picture of Enzymatic Catalysis by Hydroxynitrile Lyases from Hevea Brasiliensis: Protonation Dependent Reaction Steps and ResidueGated Movement of the Substrate and the Product. Phys. Chem. Chem. Phys. 2014, 16, 26864-26875. 36. Zhao, Y.; Chen, N.; Wang, C.; Cao, Z. A Comprehensive Understanding of Enzymatic Catalysis by Hydroxynitrile Lyases with S-Stereoselectivity from the Α/Β-Hydrolase Superfamily: Revised Role of the Active-Site Lysine and Kinetic Behavior of Substrate Delivery and Sequential Product Release. ACS Catal. 2016, 6, 2145-2157. 37. Rooklin, D. W.; Lu, M.; Zhang, Y. Revelation of a Catalytic Calcium-Binding Site Elucidates Unusual Metal Dependence of a Human Apyrase. J. Am. Chem. Soc. 2012, 134, 15595-15603. 38. Zhang, Y.; Lee, T.-S.; Yang, W. A Pseudobond Approach to Combining Quantum Mechanical and Molecular Mechanical Methods. J. Chem. Phys. 1999, 110, 46-54. 39. Chen, X.; Zhang, Y.; Zhang, J. Z. An Efficient Approach for Ab Initio Energy Calculation of Biopolymers. J. Chem. Phys. 2005, 122, 184105. 40. Zhang, Y. Pseudobond Ab Initio QM/MM Approach and Its Applications to Enzyme Reactions. Theor. Chem. Acc. 2006, 116, 43-50. 41. Zhang, Y. Improved Pseudobonds for Combined Ab Initio Quantum Mechanical/Molecular Mechanical Methods. J. Chem. Phys. 2005, 122, 024114.

ACS Paragon Plus Environment

ACS Catalysis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 18

42. Zhang, Y.; Liu, H.; Yang, W. Free Energy Calculation on Enzyme Reactions with an Efficient Iterative Procedure to Determine Minimum Energy Paths on a Combined Ab Initio QM/MM Potential Energy Surface. J. Chem. Phys. 2000, 112, 3483-3492. 43. Beeman, D. Some Multistep Methods for Use in Molecular Dynamics Calculations. J. Comput. Phys. 1976, 20, 130-139. 44. Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684-3690. 45. Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P. Advances in Methods and Algorithms in a Modern Quantum Chemistry Program Package. Phys. Chem. Chem. Phys. 2006, 8, 3172-3191. 46. Ponder, J. W. TINKER 4.2; Washington University School of Medicine, Saint Louis, MO, 2004. 47. Magistrato, A.; Sgrignani, J.; Krause, R.; Cavalli, A. Single or Multiple Access Channels to the Cyp450s Active Site? An Answer from Free Energy Simulations of the Human Aromatase Enzyme. J. Phys. Chem. Lett. 2017, 8, 2036. 48. Da, L. T.; Wang, D.; Huang, X. Dynamics of Pyrophosphate Ion Release and Its Coupled Trigger Loop Motion from Closed to Open State in Rna Polymerase Ii. J. Am. Chem. Soc. 2012, 134, 2399-2406. 49. Wang, B.; Sexton, R. E.; Feig, M. Kinetics of Nucleotide Entry into Rna Polymerase Active Site Provides Mechanism for Efficiency and Fidelity. Biochim. Biophys. Acta 2017, 1860, 482490. 50. Wu, J. W.; Wang, Z. X. The Complete Pathway for Erk2. J. Biol. Chem. 2007, 282, 27678-27684. 51. Turjanski, A. G.; Hummer, G.; Gutkind, J. S. How MitogenActivated Protein Kinases Recognize and Phosphorylate Their Targets: A QM/MM Study. J. Am. Chem. Soc. 2009, 131, 6141-6148.

Conformational changes of NahK upon substrate transportation, catalytic mechanism of GlcNAc phosphorylation, and water-assistant GlcNAc-1-P cleavage from Mg2+

ACS Paragon Plus Environment