Catalytic Mechanism of Peptidoglycan Deacetylase: A Computational

Dec 27, 2016 - Bacterial peptidoglycan deacetylase enzymes are potentially important targets for the design of new drugs. In pathogenic bacteria, they...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCB

Catalytic Mechanism of Peptidoglycan Deacetylase: A Computational Study Nicholus Bhattacharjee,†,∥ Mikolaj Feliks,†,∥ Md Munan Shaik,‡,§ and Martin J. Field*,† †

Dynamo Team/DYNAMOP Group, UMR5075, Université Grenoble I, CEA, CNRS, Institut de Biologie Structurale, 71 Avenue des Martyrs, CS 10090, 38044 Grenoble Cedex 9, France ‡ Division of Molecular Medicine, Boston Children’s Hospital, Boston, Massachusetts 02115, United States § Department of Pediatrics, Harvard Medical School, 3 Blackfan Street, Boston, Massachusetts 02115, United States ABSTRACT: Bacterial peptidoglycan deacetylase enzymes are potentially important targets for the design of new drugs. In pathogenic bacteria, they modify cell-wall peptidoglycan by removing the acetyl group, which makes the bacteria more resistant to the host’s immune response and other forms of attack, such as degradation by lysozyme. In this study, we have investigated the mechanism of reaction of acetyl removal from a model substrate, the N-acetylglucosamine/N-acetylmuramic acid dimer, by peptidogylcan deacetylase from Helicobacter pylori. For this, we employed a range of computational approaches, including molecular docking, Poisson−Boltzmann electrostatic pKa calculations, molecular dynamics simulations, and hybrid quantum chemical/molecular mechanical potential calculations, in conjunction with reaction-path-finding algorithms. The active site of this enzyme is in a region of highly negative electrostatic potential and contains a zinc dication with a bound water molecule. In the docked enzyme−substrate complex, our pKa calculations indicate that in the most stable protonation states of the active site the zinc-bound water molecule is in its hydroxide form and that the adjacent histidine residue, His247, is doubly protonated. In addition, there are one or two excess protons, with the neighboring aspartate residues, Asp12 and/or Asp199, being protonated. Overall, we find five classes of feasible reaction mechanisms, with the favored mechanism depending heavily on the protonation state of the active site. In the major one-excess-proton form, the mechanism with the lowest barrier (84 kJ mol−1) involves an initial protonation of the substrate nitrogen, followed by nucleophilic attack of the zinc-bound hydroxide and rupture of the substrate’s carbon−nitrogen bond. However, in the minor two-excess-proton form, four mechanisms are almost equienergetic (83−86 kJ mol−1), comprising both those that start with nitrogen protonation and those in which nucleophilic attack by hydroxide occurs first.



escape recognition by a host.10,11 As a result, PgdAs are of interest as targets for the design of new drugs that enhance host recognition and block immune leakage. Enzymatic deacetylation of cell-wall NAG by PgdA was first reported in Bacillus cereus.12 Since then, PgdAs have been identified and characterized from a wide range of Gram-positive as well as Gram-negative bacteria and zoonotic pathogens.13 Structural studies of PgdAs from a number of different bacteria have been reported. They indicate that the PgdAs are all metalloenzymes, within which the metal ion is coordinated by a conserved triad consisting of one aspartate and two histidine residues. This motif is common to all family 4 carbohydrate esterases (CE-4). Metal ions are indispensible for the enzymatic activity of PgdAs, as chelating agents abolish their activity.11 Maximum activity has been experimentally shown to occur with cobalt,11 but the bioavailability of different divalent cations, the

INTRODUCTION Peptidoglycan deacetylase is a member of the polysaccharide deacetylase family (pfam01522) of enzymes. It shares a conserved polysaccharide domain, also known as the NodB homology domain, with a number of other enzymes, including chitin deacetylases, acetylxylan esterases, xylanases, and rhizobial NodB chitooligosaccharide deacetylases.1,2 Peptidoglycan N-acetylglucosamine deacetylase (PgdA) (EC 3.5.1.-) catalyzes the removal of the acetyl group from the C2 atom of N-acetylglucosamine (NAG), which is a constituent of the peptidoglycan found in the cell walls of many bacteria. This modification of peptidoglycan is significant because the intracellular immune receptor nucleotide-binding oligomerization domain proteins, NOD1 and NOD2, of human cells cannot recognize modified peptidoglycan and so are unable to initiate an immune response against pathogens by either inducing autophagy or activating the NF-kB pathway.3−6 The modified peptidoglycan is also resistant to lysozyme action, thereby further contributing to pathogen survival.7−9 PgdAs are classified as virulence factors because pathogenic bacteria use them to deacetylate their own peptidoglycan as a mechanism to © XXXX American Chemical Society

Received: October 21, 2016 Revised: December 12, 2016

A

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

a previously determined lower-resolution structure of the same protein was also used (PDB code 3QBU; resolution 2.57 Å). Setup of the Continuum Electrostatic Model. Electrostatic calculations were performed in pDynamo19,20 coupled to the external Poisson−Boltzmann solver, extended-MEAD.21,22 The initial setup of the electrostatic model was performed in Chemistry at HARvard Macromolecular Mechanics (CHARMM).23,24 MM parameters, such as atomic charges and radii, were taken from the CHARMM force field.25 Hydrogen atoms were added to the crystal structure using the HBUILD routine in CHARMM. Titratable residues were protonated according to their standard protonation states at pH 7. The coordinates of the added hydrogens were geometryoptimized in CHARMM, whereas those of the heavy atoms were kept fixed. During the subsequent Poisson−Boltzmann calculations in pDynamo, the solvent was assigned an ionic strength (I) of 100 mM and a temperature (T) of 300 K. To calculate the molecular surface of the protein, a solvent probe with a radius of 1.4 Å and an ion-exclusion layer of 2.0 Å was used. The electrostatic potential was calculated using four focusing steps at resolutions of 2.0, 1.0, 0.5, and 0.25 Å. Each step employed a grid of 1213 nodes. The Metropolis Monte Carlo method, as implemented in pDynamo,26 was used to calculate the protonation states of the titratable groups. His86, His90, and Asp14 were excluded from the list of titratable groups, as these residues are involved in the coordination of the active-site zinc ion. Each MC run consisted of 500 equilibration scans, followed by 20 000 production scans. Double moves27 were applied for the pairs of groups whose absolute electrostatic interaction energy was calculated to be ≥8.4 kJ mol−1. Titration curves were plotted for a pH range of 0−14, at a resolution of 0.5 pH units. In some of our calculations, we also included the zinc-bound water explicitly and allowed it to be protonatable. In contrast to standard protein residues, however, there is uncertainty as to the reference pKaq a values of such waters, although those for the zinc-bound water in carbonic anhydrase have been estimated as ∼7.28,29 To circumvent this problem, we performed the calculations with a range of pKaq a values for this group, namely, 6−10. Setup of the Model for Molecular Docking. Molecular docking was performed in Autodock 4.1.30,31 Two separate docking simulations were done for N-methylacetamide and the NAG−NAM dimer. In the case of NAG−NAM, the anionic form of the molecule, with a deprotonated terminal carboxyl group, was taken. For each ligand, the complete tetramer of the PgdA enzyme was used as a receptor, although docking was performed only for chain D. However, during the subsequent reaction-path calculations, only the monomer with the NAG− NAM dimer docked in its active site was employed. AutoDockTools431 and custom Python and shell scripts were employed to prepare the structures of the receptor and ligands. During the docking procedure, standard parameters from the Amber force field were used to describe the receptor and ligand.32 Gasteiger charges33 were assigned to the ligand atoms. Nonpolar hydrogens were merged with their corresponding heavy atoms, which is a common strategy for docking in Autodock. Grids describing electrostatic and other nonbonding interactions were centered at the zinc ion. Each grid consisted of 40 × 40 × 40 nodes and was assigned a resolution of 0.375 Å. The grids were precalculated on Autogrid431 before docking. For each of the two ligands, a total of 100 docking runs were performed. The conformational space for each ligand was

intrinsic affinity of the enzyme, and quantum chemical (QC) calculations indicate that the preferred metal is zinc.13 Several PgdA crystal structures exist, but they are all, bar one, of the apo form of the enzyme. The exception is the structure of Bacillus subtilis PgdA (PDB ID: 1W1A), which is bound to N-acetyl-glucosamine but without the divalent ion in the active site.14 Because of this dearth of structural data, knowledge of the exact substrate specificities of PgdA and its reaction mechanism is limited, although proposals based on additional lines of evidence have been made. Thus, other zinc-dependent de-N-acetylases are known to have a basic mechanism, in which a conserved carboxylic acid residue near the zinc acts as a base to activate a water molecule, which then acts as a nucleophile to attack the acetate, generating an oxyanion tetrahedral intermediate. Similarly, mutagenesis studies of Streptococcus pneumoniae peptidoglycan GlcNAc deacetylase (SpPgdA) and in silico docking of the GlcNAc3 trimer in its active site suggest that an aspartic acid residue is the catalytic base and that the acetyl group of the substrate interacts directly with the zinc ion and a backbone nitrogen of an adjacent tyrosine residue during the reaction.11 Related mechanisms have also been proposed for two distinct acetylxylan esterases from Streptomyces lividans and Clostridium thermocellum, even though they have a preference for cobalt,15 and also for UDP-(3-O-(R-3-hydroxymyristoyl))-N-acetylglucosamine deacetylase (LpxC), which is a zinc-dependent deacetylase responsible for the biosynthesis of bacterial exopolysaccharide.16 H. pylori is a well-known risk factor for gastric carcinogenesis. Its PgdA (EC 3.5.1.10417) has been identified, characterized, and reported to contribute to the manner in which the bacterium avoids recognition by the host’s immune system.7,8,18 The high-resolution crystal structure of this PgdA shows an active site that is very similar to that of other characterized PgdAs. It contains a four-coordinate zinc ion, which is ligated by one aspartate (Asp14) and two histidine (His86 and His90) residues and a water molecule.13,18 The latter is hydrogenbonded with the residues and His247 and Asp12. On the basis of the proposals outlined above for SpPgdA, it is possible that Asp12 acts as a base that activates the water molecule, which, in turn, nucleophilically attacks the acetate of the substrate, leading to the formation of an oxyanion tetrahedral intermediate.11 In the previous experimental and computational work on the PgdA from H. pylori, we were able to locate the metals and their binding sites in the crystallographic structure18 and show that there is an intrinsic preference for zinc.13 In this article, we pursue our investigation of this enzyme by employing molecular dynamics (MD) simulations and hybrid QC/molecular mechanical (MM) calculations, in conjunction with reaction-path-finding algorithms, to study its reaction mechanism using the NAG/N-acetylmuramic acid (NAM) dimer as a substrate. To our knowledge, no kinetic data and little mechanistic data concerning the reaction exists. Our results give insights that are likely applicable to other members of the polysaccharide deacetylase family and could prove useful for the development of mechanism-based inhibitors for this critical class of metalloenzyme.



METHODS The electrostatic, docking, MD, and QC/MM models of the proteins employed in the present work were based on the crystal structure of the peptidoglycan deacetylase from H. pylori (PDB code 4LY4; resolution 2.20 Å). For comparative studies, B

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

nation states and hence with differing numbers of QC atoms and total QC charge. In the smallest QC regions, there were approximately 70 QC atoms, of which seven were link atoms, whereas the largest QC regions consisted of about 90 atoms, of which nine were link atoms. A schematic representation of the QC region is given in Figure 1.

sampled with the standard genetic algorithm implemented in Autodock 4. The receptor atoms were kept rigid, with the exception of the side chain of Trp127, which was allowed to move during the docking of the ligand. Finally, results from the docking were clustered by grouping conformations of the ligand within a root-mean-square coordinate deviation (RMSCD) of ≤1.0 Å. Relaxation of the Substrate−Enzyme Complex by MD. Because the flexibility of the receptor is somewhat limited during the docking procedure, the generated enzyme−substrate (NAG−NAM) complex was relaxed by performing a short MD simulation. During the simulation, constraints were applied to the distances and angles formed by the atoms coordinated to the zinc ion. A constraint was also applied to the distance between the oxygen atom of the water molecule coordinated to the zinc ion and the carbonyl carbon of the NAG side chain (i.e., the target of the nucleophilic attack). These constraints allowed the reaction coordinate to remain intact during the simulation while relaxing other parts of the enzyme−substrate structure. Apart from relaxing the complex, the simulation also permitted the entry of a second water molecule into the enzymatic pocket, which is vital for the QC/MM calculation. The MD model was prepared using custom Tcl scripts and visual molecular dynamics.34 The simulation was performed in NAMD.35 In the first step, the complex was solvated in a rectangular box of water molecules, using the TIP3 water model and periodic boundary conditions. The protein was separated from the boundary of the periodic cell by a 9 Å wide buffer of solvent molecules. Counterions (21 Na+) were distributed randomly around the protein to balance its nonzero charge. In the next step, the geometry of the model was optimized by using 1000 steps of the steepest descent method. The simulation was preceded by a heating phase, in which the temperature was gradually increased from 20 to 300 K, at a rate of 10 K ns−1. Both heating and relaxation were done in the NPT ensemble at a normal pressure of 1 bar and with a time step of 2 fs. The production run was performed for 5 ns. Setup of the QC/MM Model. The catalytic cycle of PgdA was studied using the hybrid QC/MM method,36−38 as implemented in pDynamo.20 The ORCA39 package, coupled to pDynamo, was employed for calculation of the density functional theory (DFT) QC energies and gradients and their associated QC/MM electrostatic interaction terms. As a starting point for the calculations, a low-energy snapshot of the enzyme−substrate complex from the previous MD simulation was selected. To reduce the computational cost, the model was truncated to a sphere of residues and water molecules within 15 Å around the center of the active site. The sphere included mostly the D chain as well as parts of the A−C chains and solvent molecules. Only complete protein residues or water molecules were incorporated into the truncated model, giving between 8033 and 8037 atoms in total (depending upon the protonation state of the protein). To keep the integrity of the model after the removal of the outer parts, the positions of the atoms at a distance of 10 Å or higher from the QC region were kept fixed. A number of QC regions for the QC/MM model were employed. The smallest consisted of the side chains of Asp12, Asp14, His86, His90, His247, the central zinc ion, two water molecules, and the reactive part of the docked substrate. In some calculations, the side chains of Asp199 and/or Asp200 were also included. As will be discussed later, we considered a number of different QC/MM systems with differing proto-

Figure 1. QC atoms employed in the QC/MM calculations of the PgdA/NAG−NAM complex. The smallest QC region contained all of the atoms shown, except for the side chains of Asp199 and Asp200. The largest QC region contained all of the atoms. Note that the protonation states of some of the residues were variable, complete details of which are given in the text.

The reaction profiles reported in the Result section were fully optimized using the largest QC/MM model and the nudged elastic band (NEB) method, as implemented in pDynamo.40,41 For these optimizations, the CHARMM27 force field was employed for the MM potential and the B3LYP DFT functional together with dispersion corrections and a triple-ζ plus polarization basis set for the QC region (ORCA keywords DFT-ENERGY D3). The final energy profiles were adjusted with a Poisson−Boltzmann solvation energy correction to estimate the neglect of the long-range interactions in the reduced QC/MM models. This was determined using the charges from the CHARMM27 force field for the MM atoms and electrostatic potential-fitted (ORCA keyword CHELPG) charges for the QC atoms, determined afresh at each point along the final NEB pathways using the full QC/MM potential. However, we note that these corrections were always small, with changes of at most ±5 kJ mol−1 with respect to the uncorrected values. The final NEB pathways were constructed using a multistep procedure. We started by geometry-optimizing the initial reactant states, in the appropriate protonation state, using structures resulting from the MD simulation. The initial optimizations were done with a pure MM potential, but these were followed by QC/MM optimizations with either a semiempirical AM1 QC method or a DFT method with the Becke−Perdew (BP) functional and a small basis set. Possible intermediates along the reaction paths were then generated from the reactant states using potential energy surface (PES) scans, in which QC/MM geometry optimizations were carried out as a function of a chosen reaction coordinate (normally a distance). These scans were performed incrementally (e.g., in steps of ±0.1 Å for a distance) until another minimum-energy C

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Comparison of the docking results for (a) N-methylacetamide and (b) the NAG−NAM dimer. For each ligand, only the three most favorable conformations are shown. The ligands and key residues in the active site are displayed in licorice and ball-and-stick representations, respectively. The catalytic zinc ion is represented as a yellow ball. Important intermolecular distances are shown as dashed lines. Small bold numbers indicate the numerical values of the distances (Å). Nonpolar hydrogens have been merged to their parent heavy atoms before the docking simulations. Both ligands bind tightly to the active site. However, binding of the NAG−NAM dimer requires rotation of the Trp127 side chain and thus opening of the active site toward the solvent so that it can accommodate bulky parts of the ligand molecule.

conformation was reached. These intermediates were themselves geometry-optimized and used as the starting points for further scans, with the procedure being repeated until the product state was reached. Once potential intermediates had been identified, partial NEB pathways were constructed from the PES scans between neighboring intermediates (including reactants and products) and were geometry-optimized. These partially optimized pathways were then joined to obtain complete pathways going from reactants to products, and these in turn were geometry-optimized using the NEB algorithm. Several NEB calculations of the complete pathways were performed using QC/MM potentials of increasing precision. Thus, initial optimizations were done with the BP functional and a small basis set for the QC region and a relatively low convergence tolerance for the NEB algorithm, followed by calculations with increasingly precise DFT methods and higher NEB tolerances until the final reported pathways were obtained. We note that during these NEB calculations the endpoint reactant and product structures were also optimized so that they too could fully adapt to the different QC/MM potentials. We also note that this scanning/partial NEB/complete NEB procedure is not guaranteed to generate feasible pathways and so many different hypotheses will often need to be tried before reasonable pathways can be found.

During the docking simulations, both ligands were found to adopt conformations suitable for catalytic reaction, as seen in Figure 2. For catalysis to occur, the oxygen atom of the zinccoordinated water should be in close contact with the carbonyl carbon atom of the substrate. In the case of N-methylacetamide, 95 of 100 docking simulations led to the ligand being closely bound to the active-site zinc ion of the protein. The side chain of Trp127 was included in the flexible region during the docking, but it consistently adopted a conformation aligning closely with that from the crystal. Only five simulations showed the ligand outside of the binding pocket. For the best-fitting conformation of N-methylacetamide in the active site of PgdA, we predicted a binding energy of −8.8 kJ mol−1, as calculated with the scoring function of Autodock 4,32 and a critical Owat··· Ccarb distance of 3.45 Å (see Figure 2). The most favorable docked conformations show the carbonyl oxygen of the ligand at a distance of approximately 4.0 Å from the zinc ion. Thus, we propose that the Zn2+···Ocarb electrostatic interaction is important for the proper orientation of the ligand in the active site. In the case of NAG−NAM, the situation becomes more complicated, as now the ligand contains a bulky group whose volume is greater than that of the binding pocket. Our initial docking simulations revealed that the binding of the dimer is possible only after the receptor becomes slightly flexible. We included the side chain of Trp127 into the flexible region on the basis of the observation that the experimentally determined flexibility of this residue, measured by its B-factor, is distinctively higher compared to that of the other residues in this region of the protein. Thus, Trp127 may play a role as a gate at the entrance to the active site. Clustering of the 100 docking simulations revealed 11 conformational clusters containing more than one member. However, the cluster with the lowest binding energy (−3.8 kJ mol−1) does not necessarily correspond to the reactive conformation of the substrate. Instead, for the three best-fitting, most reactive conformations, the binding energies were calculated to be −0.6, +3.8, and +16.7 kJ mol−1. The positive binding energy for two of these conformations suggests some strain in the enzyme−substrate complex. Hence, the geometry of the complex must be relaxed before study of the mechanism. Nonetheless, the carboxyl and



RESULTS Preparation of the Protein Models. Docking of Ligands to the Active Site. Two ligands were docked as putative substrates into the active site of PgdA, namely, Nmethylacetamide and the NAG−NAM dimer. For simplicity, we will refer to the NAG−NAM dimer simply as “the dimer”. N-Methylacetamide is the simplest molecule containing both the carboxyl and amine groups necessary for catalysis by PgdA. A docking of this ligand is instructive, as the available crystal structure of PgdA shows the protein in its closed conformation, and a small ligand, such as N-methylacetamide, may be able to fit into the rather compact active site of this form. However, the NAG−NAM dimer is a better model for the physiological substrate of the protein. D

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

MM calculations. The RMSCD of the tetramer was found to attain a stable value during course of the simulation. Protonation States of the Protein. We performed pKa calculations on the protein both before and after the docking calculations, but we have deferred discussion until now to facilitate comparison between the two. We note that the results of the pKa calculations with different pKaq a values for the zincbound water were broadly similar and also that the results were consistent between the different chains of the protein. Therefore, we focus our discussion on the calculations performed with a pKaq a value of 7 and on the D chain of the protein, unless otherwise stated. For the initial calculations on the protein without substrate, the protonatable residues in the protein were in their standard, solution states, with the exception of His54 and His147, which were doubly protonated. Of the remaining histidines, His86, His90, His151, His247, and His270 were δ-protonated and the rest were ϵ-protonated. Within the active site, the zinc is coordinated to Asp14 (deprotonated), to His86 and His90 (both δ-protonated), and to water. There was some variability in the protonation of this latter residue at the lowest pKaq a value, 6, for which we performed calculations, as the fraction of the hydroxide form increased significantly. This was accompanied by larger populations of the protonated form of Asp12 and of the doubly protonated form of His247, both of which are within the hydrogen-bonding distance of the water. Nevertheless, for the docking calculations and the subsequent MD simulations, we employed a zinc-bound water, a deprotonated Asp14, and a δ-protonated His247. The pKa calculations on the docked form of the protein resulted in identical residue protonation states to those from the previous calculations, except that now there was a clear preference for the hydroxide form of the zinc-bound water and for protonated Asp12 and doubly protonated His247. This preference persisted even in the calculations performed with the highest pKaq a value, 10, which we tried for the zinc-bound water. To illustrate possible variations in the protonations of the residues in the active site, we detail in Table 1 the lowest energy states for the D chain of PgdA docked with the NAG−NAM substrate. The proton-neutral state of the active site is one in which there are two protons (i.e., all groups in their standard, solvent reference forms). However, it can be seen that the three- and four-proton forms, with one and two excess protons, respectively, are the most stable. This can be rationalized by noting that the protein is highly negatively charged, with each monomer having a nominal, reference charge of −11 (this includes the +2 charge of the zinc cation). The most stable three- and four-proton forms both have hydroxide and protonated Asp12 (states 1 and 2, respectively). The equivalent three- and four-proton forms with water and deprotonated Asp12 (states 3 and 5, respectively) are both approximately 20 kJ mol−1 higher in energy, indicating that it is energetically more favorable to have a neutral Asp12 and a monocationic Zn−hydroxide group than an anionic Asp12 and a dicationic Zn−water group. The energy difference between the three- and four-proton forms is 9.2 kJ mol−1, indicating an approximately 97.5/2.5% population split at room temperature. The extraproton in the four-proton form occurs on Asp199 (see Table 1). Reaction-Path Calculations. We investigated the catalytic reaction for the PgdA−NAG−NAM complex in several different protonation states and found five broad categories of

amine groups of the dimer align very closely with their counterparts in the docked N-methylacetamide (see Figure 2 for a comparison of the docking of both ligands). The observation that the reacting parts occupy the same location regardless of the docked ligand may suggest that the binding pocket has a high affinity for these types of functional groups. For the best-fitting reactive conformation of NAG−NAM, the shortest Owat···Ccarb distance between the attacking water and substrate was found to be 3.29 Å. In our docking calculations and subsequent MD and reaction-path simulations, no evidence was found for the substrate significantly modifying the immediate coordination environment of the zinc ion, which includes interactions with Asp14, His86, His90, and water. Changes that could be envisaged are expansion of the coordination shell to five by the introduction of the substrate carbonyl or amine nitrogen and replacement of one of the existing groups by the substrate carbonyl. All of these possibilities appear to be either energetically or sterically disfavored. Relaxation of the Docked Complex. To further assess the binding of the ligands to the active site of PgdA, we performed an MD simulation of the NAG−NAM docked complex. Bond and angle constraints were applied between the aforementioned atoms during the 5 ns simulation to keep the zinc ion coordination intact. Apart from relaxing the enzyme−substrate complex, the main objective of the simulation was to obtain a second water molecule inside the enzymatic pocket. This is required because the crystal structure contains only one water molecule in the enzymatic pocket, bound to the zinc ion, and a second water molecule may be vital for reaction. As indicated in Figure 3, the simulation allows seepage of a second water molecule into the enzymatic pocket from the bulk solvent. This water molecule was later included in the QC region for QC/

Figure 3. Seepage of the second water molecule (red dots) toward the active site during the 5 ns constrained simulation. Atoms of the active site are shown as ball and sticks and the zinc ion as the yellow sphere. E

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 1. Energies of the First 50 Protonation States at pH 7 of the D Chain of PgdA Docked with the NAG−NAM Dimera,b state

ΔG (kJ mol−1)

water

1

0.0

−1

2

9.2

−1

3

21.2

0

4

23.4

−1

5 6 7

30.7 32.6 35.6

0 0 0

8 9

37.6 39.2

0 0

10

39.7

−1

11 12 13 14

45.6 50.4 51.9 53.6

−1 −1 −1 −1

15 16 17 18 19 20 21 22

54.4 55.5 58.9 60.1 65.0 66.3 70.4 71.7

−1 0 0 0 0 0 0 −1

23 24

73.1 74.4

0 0

25 26

75.9 79.1

0 0

state

ΔG (kJ mol−1)

water

4 3 4

27 28 29 30 31 32 33 34 35 36

79.2 80.3 81.7 81.8 83.5 87.5 93.2 94.0 94.8 95.7

0 0 −1 −1 −1 −1 −1 0 0 −1

−1 −1

3 4

37 38

103.1 103.2

0 −1

−1

−1

2

0 0 0 0

0 −1 0 −1

−1 0 0 0

3 3 4 4

39 40 41 42

105.7 105.7 106.6 110.5

−1 −1 −1 0

0 0 0 −1 −1 −1 −1 −1

0 −1 −1 −1 −1 −1 0 −1

−1 0 0 0 −1 0 0 0

3 4 4 3 2 3 4 3

43 44

126.9 132.8

−1 0

45 46 47 48 49 50

138.5 148.4 154.3 163.2 166.4 168.4

0 −1 0 0 0 −1

−1 −1

0 −1

−1 0

3 4

0 0

0 0

−1 −1

4 5

Asp199

Asp200

no. of protons

0

−1

−1

3

0

0

−1

4

−1

−1

−1

3

−1

0

−1

3

0 0 0

0 −1 −1

−1 −1 −1

−1 −1

0 0

−1

His247 Asp12 δ, ϵ, +1 δ, ϵ, +1 δ, ϵ, +1 δ, ϵ, +1 δ, 0 δ, 0 δ, ϵ, +1 δ, 0 δ, ϵ, +1 δ, ϵ, +1 δ, 0 ϵ, 0 ϵ, 0 δ, ϵ, +1 ϵ, 0 δ, 0 ϵ, 0 ϵ, 0 δ, 0 δ, 0 ϵ, 0 δ, ϵ, +1 ϵ, 0 δ, ϵ, +1 ϵ, 0 δ, ϵ, +1

His247 δ, ϵ, ϵ, δ, δ, δ, ϵ, δ, ϵ, δ,

0 0 0 0 0 0 0 0 0 ϵ, +1 ϵ, 0 δ, ϵ, +1 ϵ, 0 ϵ, 0 δ, 0 δ, ϵ, +1 δ, 0 δ, ϵ, +1 −1 δ, 0 −1 −1 −1 δ, 0

Asp12

Asp199

Asp200

no. of protons

−1 0 −1 0 0 0 0 0 0 −1

0 −1 0 −1 −1 0 −1 0 0 0

0 −1 0 −1 0 0 −1 0 0 0

4 3 3 2 3 4 2 5 5 4

−1 0

−1 0

−1 0

2 5

−1 −1 −1 0

−1 0 0 −1

0 −1 −1 0

2 2 2 5

−1 −1

0 0

0 0

3 5

0 −1 0 −1 0 −1

0 −1 −1 0 0 −1

0 0 0 0 −1 −1

4 2 3 3 3 1

a

In which the protonations of four active-site residues and the zincbound water were allowed to vary. bThe protonations of the remaining residues were kept fixed at their optimal values. The nominal charges of the residues are shown, together with the explicit protonation sites for His247. For these results, the pKaq a value of the zinc-bound water was taken to be 7.

feasible mechanisms, the schematics of which are shown in Figure 4. In each case, we started with the zinc-bound water in its hydroxide form, as the pKa calculations indicated that this is the most stable (at least when there are one to four protons in the active site). In addition, in three of the five mechanisms, the second water in the active site played a critical role. Mechanisms I−III start with nucleophilic attack of the hydroxide on the carbonyl carbon of the NAG group of the substrate, leading to a tetrahedral oxyanion. In mechanism I, this is followed by collapse of the oxyanion and internal transfer of the hydrogen from the attacking hydroxide ion to the nitrogen of NAG’s N-acetyl group. The products are an acetate ion and neutral glucosamine. Mechanism II, like mechanism I, starts with nucleophilic attack by hydroxide. In the second step, however, the nitrogen of the N-acetyl group abstracts a proton from the second water in the active site, which, in turn, deprotonates the hydroxyl group of the oxyanion species formed in the initial step. Mechanism III also involves nucleophilic attack by hydroxide, followed by a double proton transfer, but this time, the second proton comes from a neighboring protonated aspartate residue, giving acetic acid and glucosamine as products. The aspartate can be either Asp12 or Asp199, depending on the protonation state of PgdA.

In contrast to mechanisms I−III, mechanisms IV and V start with a concerted double proton transfer from a proton donor to the substrate nitrogen. The subsequent protonated substrate is then attacked nucleophilically by the zinc-bound hydroxide at the carbonyl carbon to give products. The difference between mechanisms IV and V is the initial proton donor. In mechanism IV, it is His247 and the transfer occurs via the zinc-bound hydroxide, whereas in mechanism V, it is a neighboring aspartate, either Asp12 or Asp199, and the proton transfer occurs via the second water in the active site. Altogether, we obtained over 100 well-refined pathways for these five mechanisms for the PgdA−NAG−NAM complex in various protonation states and various conformations. Although the three- and four-proton states of the active site are physiologically the most relevant ones (states 1 and 2 in Table 1), we also evaluated the pathways in the one- and twoproton states of the active site (states 10, 30, and 50 in Table 1), both for interest and to check the sensitivity of our calculations. Before further detailing the three- and four-proton pathways, we make some general remarks. First, the products were more stable than the reactants, so the overall reaction was exothermic, although the differences were greater for mechanism I and II pathways that generated an acetate ion rather than acetic acid, both of which coordinate the zinc. F

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. Schematics of the mechanisms studied for the catalytic reaction. The aspartate in mechanisms III and V can be either Asp12 or Asp199, depending upon the protonation state of PgdA.

Second, the pathways of the lowest energy had barriers of approximately 80, 80, 110, and 160 kJ mol−1 for the three-, four-, two-, and one-proton states, respectively. No mechanism was uniformly favored, and, often, several mechanisms were competitive. We note, however, that, overall, mechanism I pathways had the highest barriers, followed closely by mechanism II pathways. Third, mechanism I−III pathways always possessed two principal barriers: the first to nucleophilic attack by hydroxide, leading to a metastable oxyanion

intermediate, and the second to proton transfer. Either one of these barriers could be rate-determining, depending upon the protonation state of the complex, and the environment of the oxyanion intermediate often required small, low-barrier rearrangements for proton transfer to take place. Finally, and perhaps counterintuitively, the mechanism IV and V pathways were often as energetically favorable as those in which nucleophilic attack occurred first. G

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 5. Images of significant structures in the preferred mechanism of the most stable protonation state of the enzyme−substrate complex: (a) reactants (0 kJ mol−1); (b) transition state for the double proton transfer (84 kJ mol−1); (c) nitrogen-protonated intermediate (68 kJ mol−1); and (d) products (−18 kJ mol−1). All important residues participating in the reaction are labeled in (a), except for the zinc dication (the yellow sphere), the zinc-bound hydroxide, and the second water in the active site. Unlabeled residues that coordinate the zinc but do not participate in the reaction are Asp14 (in front of Asp12) and His86 and His90 (to the left of the zinc). Important interactions, with distances (Å), are shown as black dotted lines.

We now consider the four-proton state of the complex (state 2 in Table 1). This existed in two principal conformations, which differed depending upon whether the hydroxyl proton on the side chain of Asp199 was cis or trans to the doubly bonded oxygen. The cis conformation was favored by approximately 20 kJ mol−1, in which the Asp199 proton was strongly hydrogen-bonded to the second water molecule in the active site, whereas in the trans orientation, the proton pointed away from the active site toward residues Trp128 and Trp196. Although we determined reaction paths for the trans conformer, we found that the most favorable of these had barriers approximately 30 kJ mol−1 higher than those for the cis conformer. This, combined with the trans conformer’s higher energy of 20 kJ mol−1, means that we discuss only the cis form in what follows. In contrast to the three-proton system, the cis four-proton system exhibited four almost equienergetic mechanisms. These, in order, were mechanisms II, III (Asp199), III (Asp12), and V (Asp199), with barriers of 83, 83, 83, and 86 kJ mol−1, respectively. In the first three of these, the initial steps were identical and involved nucleophilic attack by hydroxide on the substrate. This was always rate-limiting and resulted in a tetrahedral anionic substrate transition-state structure. For mechanisms II and III (Asp199), the subsequent steps of proton transfer and carbon−nitrogen bond rupture were lower in energy and occurred simultaneously. For mechanism III

The most favorable mechanism by far for the three-proton system (state 1 in Table 1) was mechanism V (Asp12), with an overall barrier of 84 kJ mol−1. The second most favorable mechanism was mechanism IV, with a significantly higher barrier of 121 kJ mol−1. Mechanism V proceeded in two steps. The first was rate-limiting, with a barrier of 84 kJ mol−1. It involved a double proton transfer from Asp12, via the second water molecule, to the substrate nitrogen. The transition state for this step involved partially transferred protons, with significant stabilization of the second water by hydrogenbonding to Asp199. This step generated a metastable protonated substrate intermediate, with an energy 68 kJ mol−1 above that of the reactants. The second step involved simultaneous nucleophilic attack on the substrate by hydroxide and rupture of the carbon−nitrogen bond. It was an essentially barrierless process, due to the metastable nature of the intermediate, and it led to products with an energy of −18 kJ mol−1 with respect to that of the initial reactants. Images of significant structures from this pathway are given in Figure 5, which clearly show the network of interconnected hydrogen bonds between the residues that are implicated in the reaction. Perusal of the image of the product structure also suggests that a double proton transfer from the acetic acid to Asp199 via the second water molecule would be possible, thereby adjusting the protonation state of one of the deacetylation products. H

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Images of significant structures in the preferred mechanisms of the second most stable protonation state of the enzyme−substrate complex: (a) reactants (0 kJ mol−1) and (b) the tetrahedral anionic transition state that occurs in mechanisms II and III as a result of nucleophilic attack by the hydroxide on the substrate (83 kJ mol−1). See the caption of Figure 5 for nomenclature.

(Asp12), however, an intermediate was observed, 26 kJ mol−1 lower in energy than that of the transition state, followed by an 18 kJ mol−1 barrier-to-proton transfer and bond rupture. This intermediate involves a strong hydrogen bond between the second water molecule in the active site and Asp199; therefore, it is clearly more favorable to protonate the substrate from Asp199 than from Asp12. The fourth mechanism, mechanism V (Asp199), occurred in a very similar manner and with a similar barrier to that of the mechanism V pathway for the three-proton system, except that Asp199 was the proton donor. The mechanism V (Asp12) pathway for this system was also calculated and had a barrier of 115 kJ mol−1. Images of the reactant structure and the tetrahedral anionic substrate transition-state structure for mechanisms II and III are given in Figure 6 for the four-proton system. Inspection of the reactant structure helps explain why mechanisms II and III have much lower barriers than those in the three-proton system. On comparison with Figure 5a, it can be seen that the presence of the extraproton on Asp199 forces a reorganization of the hydrogen bonds of the second water. This, in turn, causes reorientation of the hydroxide, putting it in a better position for nucleophilic attack. By contrast, the extraproton does not affect the barrier for the mechanism V pathway, except for making the transfer from Asp199, instead of Asp12, more favorable. In the starting reactant structures for the mechanisms of the three- and four-proton systems (Figures 5a and 6a, respectively), the zinc is coordinated as follows: a fourcoordinate “inner” shell consisting of Asp14, His86, His90, and the hydroxide, all with zinc−ligand distances of approximately 2 Å, and a two-coordinate “outer” shell consisting of the second carboxylate oxygen from Asp14 (3.00 and 2.83 Å in Figures 5a and 6a, respectively) and the carbonyl from the substrate (3.40 and 3.03 Å in Figures 5a and 6a, respectively). Clearly, the inner-shell interaction of zinc and bound water is crucial for activating the water for reaction. Similarly, the outer-shell interaction between zinc and the substrate carbonyl helps orient the substrate in the active site and also activate the carbonyl for attack by hydroxide. To end this section, we note that we did not investigate release of products from the active site, as this is beyond the scope of our study. We envisage that the release of the glucosamine−NAM product would occur first, followed by displacement of the zinc-bound acetate or acetic acid by water. Adjustments of the protonation states of the residues in the

active site and the products, acetic acid and glucosamine, would occur in this release phase, as they become exposed to solvent.



CONCLUSIONS In this study, we have investigated the reaction mechanism of acetyl group removal from the NAG of a model substrate by the bacterial polysaccharide deacetylase from H. pylori. Because of the absence of an experimental structure of the enzyme− substrate complex, we first undertook docking calculations to determine the optimum binding modes of two model substrates, a minimal N-methylacetamide and a more realistic NAG−NAM dimer. These were combined with pKa calculations to determine the preferred protonation states of the docked enzyme−NAG−NAM complex and MD simulations to permit a proper solvation of the complex. Finally, reaction paths for a variety of possible mechanisms were determined using hybrid QC/MM potentials in conjunction with a reaction-path-finding approach. We would like to highlight a number of results arising from our study as follows: • The docking calculations identified a binding conformation for NAG−NAM that was ideal for nucleophilic attack of the zinc-bound water on the substrate’s carbonyl carbon. The validity of this structure was supported by the very similar structure found for Nmethylacetamide, even though the latter had a potentially much greater range of possible binding modes. • The enzyme has a large preponderance of negatively charged residues. As a result, the most stable protonation states of the enzyme−substrate complex were determined to be the ones in which there were either one or two excess protons in the active site, located on Asp12 and/or Asp199. Despite this, the pKa calculations showed that the most stable form of the zinc-bound water was hydroxide, with the neighboring His247 being doubly protonated and hence positively charged. • There is a well-formed and stable hydrogen-bond network that connects the principal species participating in the reactions. Thus, the zinc-bound hydroxide interacts with His247 and the second water in the active site, whereas the latter has strong interactions with Asp12 and Asp199 and, to a lesser extent, the substrate’s amide nitrogen. The reaction-path calculations indicate that proton transfer throughout this network can occur quite readily. I

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B • The preferred mechanism for the most stable protonation state of the enzyme−substrate complex is one in which protonation of the substrate nitrogen occurs first, followed closely by nucleophilic attack on the substrate by the zinc-bound hydroxide and rupture of its carbon− nitrogen bond. This result is perhaps surprising, given that it might be expected that nucleophilic attack occurs first, but can be rationalized by noting the nonoptimum orientation of the zinc-bound hydroxide for nucleophilic attack and that passage through an initial positively charged protonated substrate is not disfavored energetically, given the strong negative electrostatic potential in the active site. • Compared with the most stable protonation state, the second most stable state has an additional proton (on Asp199). This extraproton serves to reduce the negative electrostatic potential of the active site and also causes a small, but significant, reorganization of the hydrogen bonds around the second water. In particular, it reorients the zinc-bound hydroxide so that it is in a more favorable position for nucleophilic attack. As a result of these relatively subtle changes, there are now four preferred, almost equienergetic mechanisms. One of these is similar to the optimum mechanism found for the most stable protonation state of the enzyme, whereas the other three involve nucleophilic attack first, followed by protonation of the substrate nitrogen and carbon−nitrogen bond rupture. In summary, these results provide an excellent illustration of the plasticity that exists in an enzyme’s active-site structure and mechanism and of how quite a small change in the active site environment, in this case, protonation state, can switch the preferred mechanism significantly. They should also help provide insight for future studies aiming at designing inhibitors that can block the action of PgdA and similar enzymes.



(3) Meylan, E.; Tschopp, J.; Karin, M. Intracellular Pattern Recognition Receptors in the Host Response. Nature 2006, 442, 39−44. (4) Mengin-Lecreulx, D.; Lemaitre, B. Structure and Metabolism of Peptidoglycan and Molecular Requirements Allowing its Detection by the Drosophila Innate Immune System. J. Endotoxin. Res. 2005, 11, 105−111. (5) Strober, W.; Murray, P.; Kitani, A.; Watanabe, T. Signalling Pathways and Molecular Interactions of NOD1 and NOD2. Nat. Rev. Immunol. 2006, 6, 9−20. (6) Suarez, G.; Romero-Gallo, J.; Piazuelo, M.; Wang, G.; Maier, R.; Forsberg, L.; Azadi, P.; Gomez, M.; Correa, P.; Peek, R. J. Modification of Helicobacter pylori Peptidoglycan Enhances NOD1 Activation and Promotes Cancer of the Stomach. Cancer Res. 2015, 75, 1749−1759. (7) Wang, G.; Olczak, A.; Forsberg, L.; Maier, R. Oxidative Stressinduced Peptidoglycan Deacetylase in Helicobacter pylori. J. Biol Chem. 2009, 284, 6790−6800. (8) Wang, G.; Maier, S.; Lo, L.; Maier, G.; Dosi, S.; Maier, R. Peptidoglycan Deacetylation in Helicobacter pylori Contributes to Bacterial Survival by Mitigating Host Immune Responses. Infect. Immun. 2010, 78, 4660−4666. (9) Vollmer, W.; Tomasz, A. Peptidoglycan N-Acetylglucosamine Deacetylase, a Putative Virulence Factor in Streptococcus pneumoniae. Infect. Immun. 2002, 70, 7176−8. (10) Boneca, I. G.; Dussurget, O.; Cabanes, D.; Nahori, M.-A.; Sousa, S.; Lecuit, M.; Psylinakis, E.; Bouriotis, V.; Hugot, J.; Giovannini, M.; et al. A Critical Role for Peptidogly- can N-Deacetylation in Listeria Evasion from the Host Innate Immune System. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 997−1002. (11) Blair, D. E.; Schuttelkopf, A. W.; MacRae, J. I.; van Aalten, D. M. F. Structure and Metal-dependent Mechanism of Peptidoglycan Deacetylase, a Streptococcal Virulence Factor. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 15429−15434. (12) Araki, Y.; Oba, S.; Araki, S.; Ito, E. Enzymatic Deacetylation of N-Acetylglucosamine Residues in Cell Wall Peptidoglycan. J. Biochem. 1980, 88, 469−79. (13) Shaik, M. M.; Bhattacharjee, N.; Bhattacharjee, A.; Field, M. J.; Zanotti, G. Characterization of the Divalent Metal Binding Site of Bacterial Polysaccharide Deacetylase Using Crystallography and Quantum Chemical Calculations. Proteins 2014, 82, 1311−1318. (14) Blair, D. E.; van Aalten, D. Structures of Bacillus subtilis PdaA, a Family 4 Carbohydrate Esterase, and a Complex with N-Acetylglucosamine. FEBS Lett. 2004, 570, 13−9. (15) Taylor, E. J.; Gloster, T.; Turkenburg, J.; Vincent, F.; Brzozowski, A.; Dupont, C.; Shareck, F.; Centeno, M.; Prates, J.; Puchart, V.; et al. Structure and Activity of two Metal Ion-dependent Acetylxylan Esterases Involved in Plant Cell Wall Degradation Reveals a Close Similarity to Peptidoglycan Deacetylases. J. Biol. Chem. 2006, 281, 10968−75. (16) Whittington, D. A.; Rusche, K.; Shin, H.; Fierke, C.; Christianson, D. Crystal Structure of LpxC, a Zinc-dependent Deacetylase Essential for Endotoxin Biosynthesis. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 8146−50. (17) Choi, H. P.; Juarez, S.; Ciordia, S.; Fernandez, M.; Bargiela, R.; Albar, J.; Mazumdar, V.; Anton, B.; Kasif, S.; Ferrer, M.; et al. Biochemical Characterization of Hypothetical Proteins from Helicobacter pylori. PLoS One 2013, 8, No. e66605. (18) Shaik, M. M.; Cendron, L.; Percudani, R.; Zanotti, G. The Structure of Helicobacter pylori HP0310 Reveals an Atypical Peptidoglycan Deacetylase. PloS One 2011, 6, No. e19207. (19) Feliks, M.; Field, M. J. Pcetk: A pDynamo-based Toolkit for Protonation State Calculations in Proteins. J. Chem. Inf. Model. 2015, 55, 2288−2296. (20) Field, M. J. The pDynamo Program for Molecular Simulations Using Hybrid Quantum Chemical and Molecular Mechanical Potentials. J. Chem. Theory Comput. 2008, 4, 1151−1161. (21) Bashford, D. An Object-Oriented Programming Suite for Electrostatic Effects in Biological Molecules: An experience report on the MEAD project. In Scientific Computing in Object-Oriented Parallel

AUTHOR INFORMATION

Corresponding Author

*E-mail: martin.fi[email protected]. ORCID

Martin J. Field: 0000-0001-8674-7997 Author Contributions ∥

N.B. and M.F. contributed equally to this work.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS N.B., M.F., and M.J.F. thank the French National Research Agency (ANR) for support (grant numbers ANR-2011-BSV5024-04 and ANR-11-BSV5-0012). The authors also acknowledge the reviewer for helpful comments.



REFERENCES

(1) Caufrier, F.; Martinou, A.; Dupont, C.; Bouriotis, V. Carbohydrate Esterase Family 4 Enzymes: Substrate Specificity. Carbohydr. Res. 2003, 338, 687−692. (2) Kafetzopoulos, D.; Thireos, G.; Vournakis, J.; Bouriotis, V. The Primary Structure of a Fungal Chitin Deacetylase Reveals the Function for two Bacterial Gene Products. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 8005−8008. J

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Environments; Lecture Notes in Computer Science; Springer: Berlin, 1997; Vol. 1343, pp 233−240. (22) Essigke, T. A Continuum Electrostatic Approach for Calculating the Binding Energetics of Multiple Ligands. Ph.D. Thesis, University of Bayreuth, 2008. https://epub.uni-bayreuth.de/id/eprint/655. (23) 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. (24) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: the Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (25) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (26) Feliks, M. Pcetk. 2015. http://github.com/mfx9/pcetk. (27) Beroza, P.; Fredkin, D. R.; Okamura, M. Y.; Feher, G. Protonation of Interacting Residues in a Protein by a Monte Carlo Method: Application to Lysozyme and the Photosynthetic Reaction Center. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 5804−5808. (28) Silverman, D.; Lindskog, S. The Catalytic Mechanism of Carbonic Anhydrase: Implications of a Rate-limiting Protolysis of Water. Acc. Chem. Res. 1988, 21, 30−36. (29) Herr, U.; Spahl, W.; Trojandt, G.; Steglich, W.; Thaler, F.; van Eldik, R. Zinc(II) Complexes of Tripodal Peptides Mimicking the Zinc(II)-coordination Structure of Carbonic Anhydrase. Bioorg. Med. Chem. 1999, 7, 699−707. (30) Goodsell, D.; Morris, G.; Olson, A. Automated Docking of Flexible Ligands: Applications of AutoDock. J. Mol. Recognit. 1996, 9, 1−5. (31) Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J. AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 2009, 30, 2785−2791. (32) Huey, R.; Morris, G. M.; Olson, A. J.; Goodsell, D. S. A Semiempirical Free Energy Force Field with Charge-based Desolvation. J. Comput. Chem. 2007, 28, 1145−1152. (33) Gasteiger, J.; Marsili, M. Iterative Partial Equalization of Orbital Electronegativity − a Rapid Access to Atomic Charges. Tetrahedron 1980, 36, 3219−3228. (34) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (35) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (36) Field, M. J. Simulating Enzyme Reactions: Challenges and Perspectives. J. Comput. Chem. 2002, 23, 48−58. (37) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. Engl. 2009, 48, 1198−1229. (38) Acevedo, O.; Jorgensen, W. L. Advances in Quantum and Molecular Mechanical (QM/MM) Simulations for Organic and Enzymatic Reactions. Acc. Chem. Res. 2010, 43, 142−151. (39) Neese, F. The ORCA Program System. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73−78. (40) Galván, I. F.; Field, M. J. Improving the Efficiency of the NEB Reaction Path Finding Algorithm. J. Comput. Chem. 2008, 29, 139− 143. (41) Aleksandrov, A.; Field, M. A Hybrid Elastic Band String Algorithm for Studies of Enzymatic Reactions. Phys. Chem. Chem. Phys. 2012, 14, 12544−12553.

K

DOI: 10.1021/acs.jpcb.6b10625 J. Phys. Chem. B XXXX, XXX, XXX−XXX