Understanding the Catalytic Mechanism and the ... - ACS Publications

May 2, 2017 - sedolisins share a unique catalytic triad (Ser-Glu-Asp) compared with the classic ..... windows were used, and for each window, 100 ps o...
0 downloads 0 Views 3MB Size
Subscriber access provided by Eastern Michigan University | Bruce T. Halle Library

Article

Understanding the catalytic mechanism and the substrate specificity of an engineered gluten hydrolase by QM/MM MD and free energy simulations Jianzhuang Yao, Haixia Luo, and Xia Wang J. Chem. Inf. Model., Just Accepted Manuscript • Publication Date (Web): 02 May 2017 Downloaded from http://pubs.acs.org on May 2, 2017

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 free 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 accessible to all readers and 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.

Journal of Chemical Information and Modeling 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 25

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

Journal of Chemical Information and Modeling

Understanding the Catalytic Mechanism and the Substrate Specificity of an Engineered Gluten Hydrolase by QM/MM MD and Free Energy Simulations Jianzhuang Yao, a* Haixia Luo,b and Xia Wang a a

School of Biological Science and Technology, University of Jinan, Jinan 250022,

P.R. China b

Key Laboratory of Ministry of Education for Conservation and Utilization of Special

Biological Resources in the Western China, Life Science School, Ningxia University, Yinchuan 750021, P.R. China

Correspondence to: Jianzhuang Yao, Ph.D. Assistant Professor, School of Biological Science and Technology University of Jinan Jinan 250022, P.R. China Phone: +86-531-82769122 E-mail: [email protected]

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

Abstract Celiac sprue, also known as gluten-sensitive enteropathy, is a chronic disease suffered by approximately 1% population in the world. Engineered enzymes have been emerging to treat celiac disease by hydrolyzing the pathogenic peptides of gluten. For example, Kuma010 has been studied experimentally and proved to be a promising gluten hydrolase in gastric conditions. However, the detailed catalytic mechanism and the substrate specificity are still unclear. In this paper, quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) and free-energy simulations were performed to determine the catalytic mechanism and the substrate specificity and the role of the active-site residues during the reaction. The results given here demonstrate that the Kuma010 shares a similar catalytic mechanism but different substrate specificity with the wild type kumamolisin-As. Binding properties of enzyme (especially mutated residues) and substrate complex are discussed, and activation free energy barriers toward different substrates have also been examined. The computational free energy results agree reasonably with the experimental data. The strategy for developing next-generation gluten hydrolase is discussed.

ACS Paragon Plus Environment

Page 2 of 25

Page 3 of 25

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

Journal of Chemical Information and Modeling

Introduction Celiac sprue is a chronic disease suffered by approximately 1% population in the world. The basis of celiac sprue is that an immune response to gluten happens in genetically susceptible individuals, leading to a series of malnutrition-related symptoms and even lymphoma in rare cases.1 Gluten is widely prevalent in our food supply, since it is commonly found in rye, barley, and wheat. In order to treat celiac disease, kumamolisin-As was engineered as a gluten hydrolase with the seven mutations (V119D/S262K/N291D/D293T/G319S/D358G/D368H, thereafter referred as Kuma010).2, 3 As an oral-administrated pharmacological treatment, Kuma010 is aiming at detoxing the gluten by hydrolyzing proline- (Pro) and glutamine- (Gln) rich metabolic peptides in the gastric environment.2, 3 Kumamolisin-As, that belongs to the sedolisins family, is a serine carboxyl peptidase from Alicyclobacillus sendaiensis.4 Sedolisins are recently discovered serine-carboxyl peptidases that belong to the family S53 of clan SB of serine peptidases (MEROPS S53),5 and enzymes in sedolisins family share some common features.6-11 First, sedolisins are rich in conserved acidic residues (Asp and Glu); second, the maximum enzymatic activity is shown at low pH (e.g., pH3-5); third, sedolisins share a unique catalytic triad (Ser-Glu-Asp) compared to the classic serine protease (Ser-His-Asp); fourth, the general acid mechanism is adopted by sedolisins to stabilize the developing negative charge during the formation of the tetrahedral intermediate, which is significantly different from the electrostatic effect adopted by the classic serine protease. Based on the crystal structure and computational studies,4,

12

the catalytic

residues of kumamolisin-As align well and are ready for the reaction to happen. For instance, a hydrogen bond network is clearly formed in the catalytic triad (consisting of Ser278, Glu78, and Asp82). The peptide hydrolysis reactions catalyzed by kumamolisin-As consist of acylation and deacylation processes.11 In the acylation process, Glu78 in the catalytic triad of sedolisins, like the His in the catalytic triad of serine protease, acts as a general base to accept a proton from the Ser278 during the nucleophilic attack and then as a general acid to transfer a proton to the leaving group

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

(N-H) of the acylation product. In the deacylation process, the Glu78 also acts as a general base to accept a proton from a water molecule during the nucleophilic attack and then as a general acid to transfer a proton to the side chain oxygen (Oγ) atom of the catalytic Ser278 residue. According to the previous studies12, 13 with combined quantum mechanical/molecular mechianical (QM/MM) molecular dynamics (MD) and free energy simulations on the reaction of kumamolisin-As, an aspartic acid (Asp164) applies a general acid mechanism to stabilize the developing negative charge during the formation of the tetrahedral intermediate in both acylation and deacylation processes. However, questions remain as to whether the engineered kumamolisin-As uses the same catalytic mechanism as the wild type enzyme, as shown in Figure 1. Moreover, Kuma010 was specifically designed to hydrolyze the PQ motif in α-gliadin, which shows autoimmune response. As a promising gluten hydrolase, Kuma010 is expected to show broad substrate specificity toward different gliadins or PQ motif containing peptides. However, experimental kinetic data shows that Kuma010 performs better in hydrolyzing substrates containing PQL motif compared to the one with PQQ motif.2,

3

Therefore, understanding the detailed

catalytic mechanism and substrate specificity of the Kuma010 is of considerable importance in protein engineering and might be helpful for designing next generation gluten hydrolase. This computational study aimed to uncover the fundamental reaction pathway for Kuma010-catalyzed hydrolysis of peptides (PFPQPQQPF and PFPQPQLPY) by performing QM/MM MD and free energy simulations. Based on the calculations and simulations, Kuma010-catalyzed hydrolysis of peptides shares the same reaction pathway with the wild type kumamolisin-As. Moreover, the details of substrate specificity are provided through the free-energy simulations and transition state comparisons.

ACS Paragon Plus Environment

Page 4 of 25

Page 5 of 25

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

Journal of Chemical Information and Modeling

Results and discussion The catalytic mechanisms of Kuma010-catalyzed hydrolysis of both substrates (PFPQPQQPF and PFPQPQLPY) were studied computationally. The reaction pathway, including both acylation and deacylation processes, is firstly described in detail using Kuma010-PFPQPQQPF complex as an example. Then, the rate-limiting transition states for different substrates were compared to explain the substrate specificity of Kuma010. Kuma010-PFPQPQQPF binding mode. The MD simulations at the QM/MM(DFTB3:CHARMM36) level, in which the QM calculation was performed at the

third-order

self-consistent

charge

density

functional

tight-binding

(DFTB3/3ob-2-1) level14 while the MM calculation was performed using the CHARMM36 force field15, revealed that the substrate (PFPQPQQPF) occupies the active site of Kuma010 very well. The key catalytic residues of Kuma010 include the catalytic triad (consisting of Glu78, Asp82, and Ser278) and oxyanion hole (consisting of the amide groups of Ser278 and the side chain of Asp164). According to the average structure (see Figure 2A), the hydroxyl oxygen on the Ser278 side chain of Kuma010 aligns well with the carbonyl carbon on the P1-Gln6 backbone of substrate with a distance of 2.38 Å, ready for nucleophilic attack. Figure 2A also shows a hydrogen bond network was formed within the catalytic triad. Therein, Glu78 forms two strong hydrogen bonds with the side chains of Ser278 and Asp82. This observation is consistent with the results from previous studies.16-18 Meanwhile, the carbonyl oxygen of the P1-Gln6 backbone, stays in the oxyanion hole, forming hydrogen bonds with the backbone NH groups of Ser278 and side chain of Asp164. It has been suggested that the Asp164 in kumamolisin-As plays a role as general acid/base catalysts during both the acylation and deacylation processes;12, 13 a similar suggestion has also been made for Asp170 in sedolisins.16,

18

In addition to the

catalytic residues, the aromatic ring of Trp129 side chain forms a hydrophobic interaction with the substrate Pro5, as shown in Figure S5. In addition, several favorable interactions are formed between the mutated residues and substrate. For example, S73K forms a salt bridge with the C-ter of the substrate. A hydrogen bond

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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 25

network is formed among P1-Gln6, S130, and H179, which is consistent with the design strategy of Kuma010.3 In the wild type kumamolisin-As, the experimental data demonstrates that D179 is apt to stabilize a positive P1 residue (e.g. Arg).4 Therefore, a D179H mutation accounts for the switch of substrate specificity of P1 position from a positive charged residue (Arg) to a polar neutral residue (Gln). D102 may contribute electrostatic effect to stabilize the positive charge of N-ter of the substrate. In consistent with the observation by David Baker’s and his coworkers,3 D104T and D169G are to remove the steric hinder. These favorable interactions help to stabilize the peptide binding in the active site of Kuma010. Acylation reaction pathway. Starting from the Kuma010-PFPQPQQPF binding structure (Michaelis-Menten complex) as the reactant state (depicted in Figure 2A; see

Figures

S1A

and

S2A

for

atom

labels),

the

combined

QM/MM(DFTB3:CHARMM36) potential of mean force (PMF) simulations were carried out to simulate the acylation process and determine the two-dimensional (2D) free energy maps associated with two sets different reaction coordinates (as shown in Figure S3). The obtained 2D free energy maps demonstrate that the acylation process follows the classical two-step reaction pathway with a tetrahedral intermediate existing between reactant state and acyl-enzyme, well-known for peptide hydrolysis reactions catalyzed by kumomalisin-As and other carboxyl peptidases in the sedolisins family.5 The average active-site structures (i.e., TS1, TI1, TS2, and AE) at the different stages of the acylation process for Kuma010 are plotted in Figure 2B-D. The acylation process is mainly reflected by the formation of a new covalent bond (C−Oγ) between the carbonyl carbon (C) on the P1-Gln6 backbone of substrate and the hydroxyl oxygen (Oγ) on Ser278 side chain of Kuma010 (first step) and breaking of a covalent bond (C−N) between the carbonyl carbon (C) and the peptide nitrogen (N) on the P1-Gln6 backbone of substrate (second step). Meanwhile, two proton transfer reactions happen associated with these two steps. The proton transfer reactions in the first step include the hydrogen (Hγ) of Ser278 gradually moves toward the oxygen (OE) atom of the Glu78 side chain and the hydrogen (HD) of Asp164 gradually moves

ACS Paragon Plus Environment

Page 7 of 25

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

Journal of Chemical Information and Modeling

toward the oxygen (O) atom of the P1-Gln6 backbone of substrate. Thus, these two residues play the important roles of the general base and acid, respectively, during the formation of TI1, consistent with the previous results on the acylation reaction of the wild type kumamolisin-As.12 As shown in Fig. S1A, five bond lengths were used to generate two reaction coordinates, denoted as RC1 = r(Oγ...Hγ) - r(OE...Hγ) - r(C...Oγ) and RC2 = r(OD...HD) - r(O... HD), for the QM/MM 2D-PMF calculations. The proton transfer reactions in the second step include the hydrogen (Hγ) of Ser278 (which now resides at Glu78) gradually moves towards the peptide nitrogen (N) on the P1-Gln6 backbone of substrate, and the hydrogen (HD) of Glu164, which now resides at the P1-Gln6 backbone of substrate, gradually moves back to the oxygen (OD) atom of Asp164. As shown in Fig. S1A, five bond lengths were used to generate two reaction coordinates, denoted as RC3 = r(N...C) + r(OE...Hγ) - r(N...Hγ) and RC2, for the QM/MM 2D-PMF calculations. Deacylation reaction pathway. Starting from the AE state (depicted in Figure 3A;

Figures

S1A

and

S2A

for

atom

labels),

the

combined

QM/MM(DFTB3:CHARMM36) potential of mean force (PMF) simulations were carried out to simulate the deacylation process and determine 2D free energy maps associated with two sets different reaction coordinates (as shown in Figure S3C-D). The obtained 2D free energy maps demonstrate that the deacylation process also follows the well-known two-step reaction pathway known for peptide hydrolysis reactions catalyzed by kumomalisin-As and other carboxyl peptidases in the sedolisins family.5 The average active-site structure for the AE-H2O complex obtained from the QM/MM MD simulations is shown in Figure 3A. As is evident from Figure 3A, the nucleophile water is well positioned to attack the carbonyl carbon of the P1-Gln6 during the MD simulations. Glu78 and Asp164 are located at the right positions to function as the general base and acid catalysts during the nucleophilic attack, respectively, which is consistent with the previous study of wild type kumomalisin-As

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

deacylation reaction.16 Asp82 forms a strong hydrogen bond with Glu78 in the acyl-enzyme complex. Figure 3B-D shows the chemical transformation from the acyl-enzyme (AE+H2O) to a tetrahedral intermediate (TI2) through transition state TS3, the oxygen (Ow) atom of the nucleophilic water gradually approaches the carbonyl carbon (C) atom while a proton (Hw) of the water molecule gradually transfers to the OE atom of Glu78 side chain. During the next reaction step from TI2 to the product state (PC depicted in Figure 1), the proton (Hw) gradually transfers from the OE atom of Glu78 side chain to the hydroxyl oxygen (Oγ) atom of Ser278 side chain while the covalent bond C−Oγ gradually breaks. As shown in Fig. S2A, for 2D-PMF calculations, the reaction coordinates associated with TS2 were RC4 = r(Ow...Hw) - r(OE...Hw) r(C...Ow) and RC2, and the reaction coordinates associated with TS3 was RC5 = r(C...Oγ) + r(OE...Hw) - r(Oγ...Hw) and RC2. Overall free energy profiles and substrate specificity in comparison with available experimental data. Depicted in Figure 4A-B is the minimum free energy profiles plotted based on the 2D free energy maps (in Figures S3-4). According to Figure 4A, the entire reaction of Kuma010 catalyzed peptide (PFPQPQQPF) hydrolysis contains acylation and deacylation processes and follows the classical catalytic mechanism of wild type kumomalisin-As and other carboxyl peptidases in the sedolisins family.5 Both acylation and deacylation processes are two-steps reactions with a tetrahedral intermediate sandwiched by two transition states. The free energy barriers (PFPQPQQPF hydrolysis) associated with the four transition states (TS1 to 4) are 10.1, 19.3, 14.9, and 14.3 kcal/mol, respectively. Thus, the second reaction step (acylation) is rate-determining for Kuma010-catalyzed hydrolysis of the peptide (PFPQPQQPF), which is also in agreement with the wild type kumamolisin-As catalyzed reaction.12, 13, 19 To elucidate the substrate specificity of Kuma010, the same QM/MM MD and free energy (PMF) simulations were performed for the rate-limiting acylation process of Kuma010-catalyzed hydrolysis of the peptide (PFPQPQLPY). Based on our simulation results, the reaction pathways of acylation processes for both substrates catalyzed by Kuma010 are identical, with the second TS

ACS Paragon Plus Environment

Page 8 of 25

Page 9 of 25

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

Journal of Chemical Information and Modeling

as the rate-limiting step. According to Figure 4B, the activation free energy barriers of PFPQPQLPY hydrolysis is 19.0 kcal/mol, which is slightly smaller than PFPQPQQPF hydrolysis catalyzed by Kuma010. According to the Figure S6, the error bars are quite small for both rate-limiting steps (PFPQPQQPF: 19.32±0.06 kcal/mol; PFPQPQLPY: 19.01±0.01 kcal/mol), showing that the calculated activation free energy barrier difference (0.3 kcal/mol) is statistically meaningful. The kcat values are not directly available for both substrates, but can be calculated based on the Michaelis-Menten curves obtained by experimental estimation2. As shown in the Fig. S1 of the Reference,2 the Vmax for PFPQPQLPY and PFPQPQQPF are roughly 3.0×10-7 M/S and 0.7×10-7 M/S, and the concentration of Kuma010 are 20 nM for both substrates. Thus, the kcat values for PFPQPQLPY and PFPQPQQPF are around 15.0 S-1 and 3.5 S-1. According to the conventional transition state theory,20 the activation free energy barriers (under 37℃) for PFPQPQLPY and PFPQPQQPF are 16.5 kcal/mol and 17.4 kcal/mol, respectively. As shown in Figure S7, DFTB3/MM tends to overestimate the rate-limiting potential energy barrier, compared to B3LYP/MM method. This explains the reason why our computational activation free energies are systematically overestimated compared to experimental activation free energy barriers. Therefore, both computational and experimental derived activation free energy barriers reasonably agree with each other, which suggests that the obtained mechanistic insights are reasonable. Depicted in Figure 4C-D are average structures of the approximate transition states obtained from the PMF calculations. Both structures show that Kuma010 provides two hydrophobic residues (L81 and I274) to stabilize P1’ residues (L and Q) of the substrates (PFPQPQLPY and PFPQPQQPF). In Figure 4D, P1’-Leu7 of the substrate (PFPQPQLPY) is stabilized by a hydrophobic cave, consisting of Leu81, Ile274, and hydrophobic moiety of P3’-Tyr9 side chain. The shortest average distances between P1’-Leu7 and Leu81, Ile274, and hydrophobic moiety of P3’-Tyr9 side chain are 3.72 Å, 3.98 Å, and 3.90 Å, respectively. In Figure 4C, the P1’-Gln7 of the substrate (PFPQPQQPF) also locates in the same hydrophobic cave consisting of Leu81 and Ile274, for the reaction to happen. The shortest distances between the

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

P1’-Gln7 and Leu81 and Ile274 are 3.38 Å and 3.71 Å. These close contacts between hydrophobic residues and hydrophilic residues is not favorable to stabilize the substrate (PFPQPQQPF) in the active site. Interestingly, the P3’-Phe9 does not directly contact with P1’-Gln7 (with a distance of 4.50 Å), which is different with PFPQPQLPY transition state complex. Instead, a water molecule moves into the active site and forms a hydrogen bond (hydrophilic interaction) with P1’-Gln7 side chain. Based on the transition states comparison, Kuma010 seems to be able to provide more favorable hydrophobic interaction to stabilize the P1’ residue (P1’-Leu7) of PFPQPQLPY rather than PFPQPQQLF (P1’-Gln7). Thus, our computational transition state structure implicates that the catalytic efficiency of Kuma010-catalyzed hydrolysis of PFPQPQQLF can be further improved by introducing polar interaction between the active site residues and the P1 residue (Gln) of the substrate. Indeed, another gluten hydrolase (named as Kuma030) with a mutation (I274T) designed by the same research group2 shows improved catalytic efficiency toward PFPQPQQLF. Benchmark calculations. The application of SCC-DFTB to describe the reactions catalyzed by sedolisins family enzymes has been proved by previous studies.12, 13, 16-19, 21, 22 Here, a benchmark calculation23, 24 was applied to compare QM/MM(DFTB3:CHARMM36)

with

QM/MM(B3LYP/6-31G*:CHARMM36)

method to prove that the DFTB3 method is valid on simulation of the peptide (PFPQPQQPF) hydrolysis reactions catalyzed by Kuma010. The results of the reaction processes comparison are depicted in Figure S1-2 and S7. The initial structures for the optimizations were obtained from the last snapshots of tetrahedral intermediates (TI1 and TI2 for acylation and deacylatoin, respectively) from the QM/MM(DFTB3:CHARMM36) PMF simulations. The optimized structure of the reactant complex from QM/MM(DFTB3:CHARMM36) shows generally good agreement with that obtained from QM/MM(B3LYP/6-31G*:CHARMM36). The hydrogen bond networks of the catalytic triad (consisting of Glu78, Asp82, and Ser278) and the oxyanion hole (Asp164 and amide group of Ser278) are aligned very well for the both QM/MM methods. The DFTB3/MM and B3LYP/MM results show the similar trends and catalytic mechanism for both acylation and deacylation

ACS Paragon Plus Environment

Page 10 of 25

Page 11 of 25

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

Journal of Chemical Information and Modeling

processes from the adiabatic mapping calculations as shown in Figure S1-2 and S7. For example, the developing negative charge of carbonyl O along with the reaction is mainly stabilized by Asp164 through a general acid/base mechanism. Although the distances optimized by two QM/MM methods are slightly different and DFTB3/MM tends to overestimate the potential energy barrier compared to B3LYP/MM, the comparison indicates that overall the DFTB3/MM Hamiltonian is able to provide a reliable description for the structure and potential energy features, and reaction pathways of Kuma010 catalyzed hydrolysis of its substrate (PFPQPQQPF).

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

Conclusion All of the QM/MM reaction-coordinate calculations, MD and PMF simulations have consistently demonstrated that the reaction pathway of Kuma010-catalyzed hydrolysis of peptides follows the classical catalytic mechanism of the serine carboxyl peptidases, consisting of the acylation and the deacylation processes. Both acylation and deacylation processes are two steps reactions with tetrahedral intermediates between reaction state and acyl-enzyme (acylation), and between acyl-enzyme and the product state (deacylation). The second transition state (TS2) in the acylation process is rate-determining for the enzymatic hydrolysis of peptides. The hydrophobic residues (e.g., Leu81 and Ile274) are suggested to account for the substrate specificity toward two substrates (PFPQPQQPF and PFPQPQLPY). This computational activation free energy barriers reasonably agree with the experiment study2, suggesting that the obtained mechanistic insights are reasonable.

ACS Paragon Plus Environment

Page 12 of 25

Page 13 of 25

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

Journal of Chemical Information and Modeling

Methods Construction of the initial Kuma010-substrates binding structures. The binding model of Kuma010 and substrates were built using a similar Kuma010-design protocol described previously by Gordon et al.3 It must be a safe and quick to generate a working model of the Kuma010 reactant complex by following the original design protocol. The crystal structure of apo-form Kuma010 (PDB ID 4NE7)2 was not used as the starting structure, because the side chain of catalytic Ser278 adopts a completely different orientation, without forming the hydrogen bond network with other residues in the catalytic traid. The crystal structure of pro-Kumamolisin (PDB ID: 1T1E)10 was used to identify the binding site of the substrates (PFPQPQQPF and PFPQPQLPY). The corresponding residues in 1T1E binding site were superimposed into the apoprotein of kumamolisin-As (PDB ID: 1SIO).4 The kumamolisin-As and substrate residues were manually mutated to the Kuma010-PFPQPQQPF and Kuma010-PFPQPQLPY complexes. The residue number of the Kuma010 (V119D/S262K/N291D/D293T/G319S/D358G/D368H)

is

based

on

the

pro-kumamolisin-As. The V119D, which locates in the propeptide domain, will be cut off during the self-activation process of Kuma010, and thus only six mutations are kept in mature Kuma010.3 According to the traditional residue number of the mature kumamolisin-As,

the

six

mutations

should

be

S73K/N102D/D104T/G130S/D169G/D179H. Hydrogen atoms of Kuma010 were added by the HBUILD module25 implemented in CHARMM.26 All acidic and basic amino acids protonation states were determined by surrounding environment under pH 4.0 condition. All protonation states are confirmed by H++.27 The Kuma010-substrates complexes were solvated by a water droplet with 22 Å radius by the standard superimposing protocol at the center of the oxygen atom (Oγ) of Ser278 side chain, and solvent water molecules within 2.8 Å of any crystal atoms were removed. A modified TIP3P water model15, 28-30 was employed for the solvent. Included in the QM-treated region of reaction processes were the side chains of Glu78, Asp82, Asp164, and Ser287, as well as a part of the substrate peptide containing the carbonyl of Pro5, the backbone of Gln6, and the Cα

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

and amide group of Gln7 or Leu7 depending on the substrates. The remaining part of the system was included in the MM-treated region. To simulate deacylation process, acylation product (QPF or LPY) was removed from the system and a water molecule was added into the QM-treated region used in the simulation of the acylation process. The divided frontier charge (DIV) link-atom scheme26, 31, 32 was applied to the QM and MM boundary treatment. The DFTB3 method implemented in the CHARMM33 was used for the QM atoms and the all-hydrogen CHARMM potential function (PARAM36)15, 28 was used for the MM atoms. The non-bonded interaction truncation cut-off was 13 Å. Stochastic boundary molecular dynamics (MD) method34 implemented in CHARMM was used with the oxygen atom (Oγ) of Ser278 side chain as the reference centre. The reaction region was a sphere with a radius (r) of 20 Å, and the buffer region extended over 20 Å ≤ r ≤ 22 Å. Newtonian equations-of-motion were solved for the reaction region (within 20 Å), and Langevin equations-of-motion were solved for the buffer regions (20−22 Å) with a temperature bath of 300 K. The temperature was controlled by Langevin thermostat.35 All atoms located longer than 22 Å away from the hydroxyl oxygen on Ser278 side chain of Kuma010 were fixed. All bonds involving hydrogen atoms were constrained by using the SHAKE algorithm.36 The initial structure for the entire stochastic boundary system was optimized using the steepest descent (SD) and adopted-basis Newton-Raphson (ABNR) methods. A 1-fs time step was used for integration of the equation of motion. The system was gradually heated from 50.0 K to 298.15 K in 100 ps and then 500 ps production run was performed for reactant complexes (RC) of Kuma010- PFPQPQQPF and Kuma010- PFPQPQLPY. Potential mean force (PMF) free energy simulation. The PMF free energy simulation was based on the QM/MM(DFTB3:CHARMM36) MD simulations. The umbrella sampling method37 implemented in the CHARMM program along with the Weighted Histogram Analysis Method (WHAM)38 was used to determine the change of the free energy (PMF) surface as a function of the reaction coordinates. The potential energy maps were first generated for both the acylation and deacylation processes by adiabatic mapping calculations starting from the last snapshot of the

ACS Paragon Plus Environment

Page 14 of 25

Page 15 of 25

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

Journal of Chemical Information and Modeling

QM/MM(DFTB3:CHARMM36) MD simulation on the reaction system. For the free energy maps (including the acylation and deacylation processes), a total of more than three thousand windows were used and, for each window, 100 ps MD simulation was performed with 50 ps for equilibration, and 90ps MD simulation with 50 ps for equilibration was applied for error bar analysis purpose. One snapshot was save every 500fs, so 100 snapshots are saved for one windows. More than 2600 windows of MD simulations were performed. The force constant of the harmonic biasing potential used in all of the PMF simulations was 150 kcal mol-1 Å-2. Benchmark calculations. The final TI1 and TI2 structures mentioned above was used to perform QM/MM reaction-coordinate calculations at the QM/MM(DFTB3:CHARMM36) level in which the QM calculation was performed at the DFTB3 level14 while the MM calculation was performed using the CHARMM36 force field15, 28 implemented in the CHARMM program. Further, for validation of DFTB3, the QM/MM reaction-coordinate calculations were carried out using B3LYP/6-31G* instead of DFTB3. The QM/MM reaction-coordinate calculations were carried out by using the GAMESS-US program39 interfaced with the CHARMM program.26, 40 Reaction coordinates used in the QM/MM reaction-coordinate calculations were defined as linear combinations of key internuclear distances involved in the reaction process. The potential energy surface was generated by the adiabatic mapping calculations starting from the last snapshot of QM/MM(DFTB3:CHARMM36) PMF simulation on the TI1 for acylation and TI2 for deacylation. The adiabatic mapping calculations at all QM/MM levels were carried out using a uniformed protocol. Specifically, a harmonic restraint with a force constant of 10000 kcal mol-1 Å-2 was added to the reaction coordinate to guide the Kuma010 acylation and deacylation reaction processes. The reaction coordinate value was increased or decreased stepwise along the reaction path (from TI1 to RC and AE for acylation or from TI2 to AE+H2O and PC for deacylation), with a step size of 0.1 Å. An adopted basis Newton Raphson (ABNR) minimization was carried out under the restraint. The reaction coordinate value then decreased or increased from RC and AE to TI1 for acylation or from

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

AE+H2O and PC to TI2 for deacylation. The forward and backward processes41 were repeated until a smooth and more reliable potential energy surface profile along the reaction path was obtained.

Supporting Information Additional Figures (S1 to 7). This information is available free of charge via the Internet at http://pubs.acs.org. Corresponding Author *Tel.: +86-513-82769122. E-mail: [email protected]. ORCID Jianzhuang Yao: 0000-0002-5017-3677 Funding This work was supported in part by the Ningxia Natural Science Foundation Project (NZ14032) and Ningxia Hui Autonomous Region Overseas Students Innovation Merit Founded project Ning Ren social letter (2014-486) and Initial Scientific Research Found of Ningxia University Talent Introduction to HXL. Notes The authors declare no competing financial interest.

ACS Paragon Plus Environment

Page 16 of 25

Page 17 of 25

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

Journal of Chemical Information and Modeling

Reference 1. Shan, L.; Molberg, Ø.; Parrot, I.; Hausch, F.; Filiz, F.; Gray, G. M.; Sollid, L. M.; Khosla, C., Structural Basis for Gluten Intolerance in Celiac Sprue. Science 2002, 297, 2275-2279. 2. Wolf, C.; Siegel, J. B.; Tinberg, C.; Camarca, A.; Gianfrani, C.; Paski, S.; Guan, R.; Montelione, G.; Baker, D.; Pultz, I. S., Engineering of Kuma030: A Gliadin Peptidase That Rapidly Degrades Immunogenic Gliadin Peptides in Gastric Conditions. J. Am. Chem. Soc. 2015, 137, 13106-13113. 3. Gordon, S. R.; Stanley, E. J.; Wolf, S.; Toland, A.; Wu, S. J.; Hadidi, D.; Mills, J. H.; Baker, D.; Pultz, I. S.; Siegel, J. B., Computational Design of an α-Gliadin Peptidase. J. Am. Chem. Soc. 2012, 134, 20513-20520. 4. Wlodawer, A.; Li, M.; Gustchina, A.; Tsuruoka, N.; Ashida, M.; Minakata, H.; Oyama, H.; Oda, K.; Nishino, T.; Nakayama, T., Crystallographic and Biochemical Investigations of Kumamolisin-As, a Serine-Carboxyl Peptidase with Collagenase Activity. J. Biol. Chem. 2004, 279, 21500-21510. 5. Oda, K., New Families of Carboxyl Peptidases: Serine-Carboxyl Peptidases and Glutamic Peptidases. J. Biochem. 2012, 151, 13-25. 6. Takahashi, K.; Oda, K., Structural Evidence That Scytalidolisin (Formerly Scytalidopepsin A) is a Serine-Carboxyl Peptidase of the Sedolisin Family. Biosci., Biotechnol., Biochem. 2008, 72, 2239-42. 7. Siezen, R. J.; Renckens, B.; Boekhorst, J., Evolution of Prokaryotic Subtilases: Genome-Wide Analysis Reveals Novel Subfamilies with Different Catalytic Residues. Proteins 2007, 67, 681-94. 8. Reichard, U.; Lechenne, B.; Asif, A. R.; Streit, F.; Grouzmann, E.; Jousson, O.; Monod, M., Sedolisins, a New Class of Secreted Proteases from Aspergillus Fumigatus with Endoprotease or Tripeptidyl-Peptidase Activity at Acidic PHs. Appl. Environ. Microbiol. 2006, 72, 1739-48. 9. Suzuki, N.; Nishibori, K.; Oodaira, Y.; Kitamura, S.; Michigami, K.; Nagata, K.; Tatara, Y.; Lee, B. R.; Ichishima, E., Grifolisin, a Member of the Sedolisin Family Produced by the Fungus Grifola Frondosa. Phytochemistry 2005, 66, 983-90. 10. Comellas-Bigler, M.; Maskos, K.; Huber, R.; Oyama, H.; Oda, K.; Bode, W., 1.2 A Crystal Structure of the Serine Carboxyl Proteinase Pro-Kumamolisin; Structure of an Intact Pro-Subtilase. Structure 2004, 12, 1313-23. 11. Wlodawer, A.; Li, M.; Gustchina, A.; Oyama, H.; Dunn, B. M.; Oda, K., Structural and enzymatic properties of the sedolisin family of serine-carboxyl peptidases. Acta Biochim Pol 2003, 50, 81-102. 12. Xu, Q.; Guo, H. B.; Wlodawer, A.; Nakayama, T.; Guo, H., The QM/MM Molecular Dynamics and Free Energy Simulations of the Acylation Reaction Catalyzed by the Serine-Carboxyl Peptidase Kumamolisin-As. Biochemistry 2007, 46, 3784-3792. 13. Xu, Q.; Li, L.; Guo, H., Understanding the Mechanism of Deacylation Reaction Catalyzed by the Serine Carboxyl Peptidase Kumamolisin-As: Insights from QM/MM Free Energy Simulations. J. Phys. Chem. B, 2010, 114, 10594–10600. 14. Gaus, M.; Cui, Q.; Elstner, M., DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931-948.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

15. Huang, J.; MacKerell, A. D., CHARMM36 All-Atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135-2145. 16. Yao, J.; Xu, Q.; Guo, H., QM/MM and Free-Energy Simulations of Deacylation Reaction Catalysed by Sedolisin, a Serine-Carboxyl Peptidase. Mol. Simul. 2013, 39, 206-213. 17. Yao, J.; Wlodawer, A.; Guo, H., Understanding the Autocatalytic Process of Pro-Kumamolisin Activation from Molecular Dynamics and Quantum Mechanical/Molecular Mechanical (QM/MM) Free-Energy Simulations. Chem. - Eur. J. 2013, 19, 10849-10852. 18. Xu, Q.; Yao, J.; Wlodawer, A.; Guo, H., Clarification of the Mechanism of Acylation Reaction and Origin of Substrate Specificity of the Serine-Carboxyl Peptidase Sedolisin through QM/MM Free Energy Simulations. J. Phys. Chem. B 2011, 115, 2470-2476. 19. Xu, Q.; Guo, H.; Wlodawer, A., The Importance of Dynamics in Substrate-Assisted Catalysis and Specificity. J. Am. Chem. Soc. 2006, 128, 5994-5995. 20. Garrett, B. C.; Truhlar, D. G., Transition State Theory in Encyclopedia of Computational Chemistry. John Wiley & Sons: Chichester, UK, 1998. 21. Guo, H. B.; Wlodawer, A.; Nakayama, T.; Xu, Q.; Guo, H., Catalytic Role of Proton Transfers in the Formation of a Tetrahedral Adduct in a Serine Carboxyl Peptidase. Biochemistry 2006, 45, 9129-9137. 22. Guo, H. B.; Wlodawer, A.; Guo, H., A General Acid-Base Mechanism for the Stabilization of a Tetrahedral Adduct in a Serine-Carboxyl Peptidase: A Computational Study. J. Am. Chem. Soc. 2005, 127, 15662-15663. 23. Yao, J.; Guo, H.; Chaiprasongsuk, M.; Zhao, N.; Chen, F.; Yang, X.; Guo, H., Substrate-Assisted Catalysis in the Reaction Catalyzed by Salicylic Acid Binding Protein 2 (SABP2), a Potential Mechanism of Substrate Discrimination for Some Promiscuous Enzymes. Biochemistry 2015, 54, 5366-5375. 24. Hou, G.; Cui, Q., Stabilization of Different Types of Transition States in a Single Enzyme Active Site: QM/MM Analysis of Enzymes in the Alkaline Phosphatase Superfamily. J. Am. Chem. Soc. 2013, 135, 10457-10469. 25. Brünger, A. T.; Karplus, M., Polar Hydrogen Positions in Proteins: Empirical Energy Placement and Neutron Diffraction Comparison. Proteins: Struct., Funct., Bioinf. 1988, 4, 148-156. 26. Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M., CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545-1614. 27. 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. 28. MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin,

ACS Paragon Plus Environment

Page 18 of 25

Page 19 of 25

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

Journal of Chemical Information and Modeling

D.; Karplus, M., All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586-3616. 29. Jorgensen, W. L., Quantum and Statistical Mechanical Studies of Liquids .10. Transferable Intermolecular Potential Functions for Water, Alcohols, and Ethers - Application to Liquid Water. J. Am. Chem. Soc. 1981, 103, 335-340. 30. Neria, E.; Fischer, S.; Karplus, M., Simulation of Activation Free Energies in Molecular Systems. J. Chem. Phys. 1996, 105, 1902-1921. 31. Field, M. J.; Bash, P. A.; Karplus, M., A Combined Quantum-Mechanical and Molecular Mechanical Potential for Molecular-Dynamics Simulations. J. Comput. Chem. 1990, 11, 700-733. 32. König, P. H.; Hoffmann, M.; Frauenheim, T.; Cui, Q., A Critical Evaluation of Different QM/MM Frontier Treatments with SCC-DFTB as the QM Method. J. Phys. Chem. B 2005, 109, 9082-9095. 33. Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M., A QM/MM Implementation of the Self-Consistent Charge Density Functional Tight Binding (SCC-DFTB) Method. J. Phys. Chem. B 2000, 105, 569-585. 34. Brooks, C. L.; Brunger, A.; Karplus, M., Active-Site Dynamics in Protein Molecules - A Stochastic Boundary Molecular-Dynamics Approach. Biopolymers 1985, 24, 843-865. 35. Adelman, S. A.; Doll, J. D., Generalized Langevin Equation Approach For Atom-Solid-Surface Scattering - General Formulation For Classical Scattering Off Harmonic Solids. J. Chem. Phys. 1976, 64, 2375-2388. 36. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C., Numerical-Integration of Cartesian Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327-341. 37. Torrie, G. M.; Valleau, J. P., Monte-Carlo Free-Energy Estimates Using Non-Boltzmann Sampling - Application to Subcritical Lennard-Jones Fluid. Chem. Phys. Lett. 1974, 28, 578-581. 38. Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M., The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules .1. The Method. J. Comput. Chem. 1992, 13, 1011-1021. 39. Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., General Atomic and Molecular Electronic-Structure System. J. Comput. Chem. 1993, 14, 1347-1363. 40. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., CHARMM - A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187-217. 41. 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.

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

Figure 1. Catalytic Mechanism of Kuma010 with a PQQ-containing substrate (PFPQPQQPF) and the role of key active-site residues for both acylation and deacylation processes. Figure 2. The QM/MM(DFTB3:CHARMM36) MD-simulated average structures of acylation process (Kuma010-catalyzed hydrolysis of a PQQ-containing substrate PFPQPQQPF) in different states: (A) RC; (B) TS1; (C) TI1; (D) TS2. For RC, all of the collected snapshots of the last 500 ps (one snapshot for 0.5 ps) of the MD trajectory were used to obtain the average structure. For all other acylation states, all of the collected snapshots of the last 50 ps (one snapshot for 0.5 ps) of the PMF trajectory were used to obtain the average structure. The distances are given in Å. Figure 3. The QM/MM(DFTB3:CHARMM36) MD-simulated average structures of deacylation process (Kuma010-catalyzed hydrolysis of a PQQ-containing substrate PFPQPQQPF) in different states: (A) AE+H2O; (B) TS3; (C) TI2; (D) TS4. For all deacylation states, all of the collected snapshots of the last 50 ps (one snapshot for 0.5 ps) of the PMF trajectory were used to obtain the average structure. The distances are given in Å. Figure 4. Free energy profiles and transition states stabilization. Minimum free energy profiles determined by the QM/MM(DFTB3:CHARMM36) calculations based two-dimensional PMF free energy maps for (A) the entire reaction process (acylation and deacylation) of Kuma010-catalyzed hydrolysis of a PQQ-containing substrate (PFPQPQQPF), and (B) The acylation process of Kuma010-catalyzed hydrolysis of a PQL-containing substrate (PFPQPQLPY). Approximate TS2 of Kuma010-catalyzed hydrolysis of (C) a PQQ-containing substrate (PFPQPQQPF) and (D) a PQL-containing substrate (PFPQPQLPY).

ACS Paragon Plus Environment

Page 20 of 25

Page 21 of 25

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

Journal of Chemical Information and Modeling

Figure 1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

Figure 2

ACS Paragon Plus Environment

Page 22 of 25

Page 23 of 25

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

Journal of Chemical Information and Modeling

Figure 3

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

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

Figure 4

ACS Paragon Plus Environment

Page 24 of 25

Page 25 of 25

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

Journal of Chemical Information and Modeling

For Table of Contents Only

ACS Paragon Plus Environment