QM/MM analysis suggests that Alkaline Phosphatase (AP) and Nucleotide pyrophosphatase/phosphodiesterase slightly tighten the transition state for phosphate diester hydrolysis relative to solution: implication for catalytic promiscuity in the AP superfamily

A. Diester hydrolysis in solution with the SCC-DFTBPR based implicit solvent model

As an important benchmark and reference, we first study the hydrolysis of MpNPP⁻ and two of its analogs (Fig.2) in solution, which have been thoroughly studied experimentally. To focus on the quality of the QM model rather than other technical details such as QM/MM coupling and sampling, we use the implicit solvent model that we have implemented and parameterized for SCC-DFTBPR.⁴¹ In this model, the solute radii are dependent on the charge distribution, which makes it particularly useful for studying solution reactions that involve highly charged species; our previous benchmark calculations suggest that the method has comparable accuracy as the SM6 model⁴², while being much more efficient (due to the use of SCC-DFTB) and having only a small number of parameters.

The potential energy surface (PES) relevant to the hydrolysis reaction is first explored by adiabatic mapping calculations, in which the reaction coordinates are the P-O^lg and P-O^nu distances; here hydroxide is the nucleophile, and “O^lg” and “O^nu” indicate the reactive oxygen in the leaving group and nucleophile, respectively. Each point in the 2D PES is obtained by starting the constrained optimization from several different initial structures and taking the lowest energy value. Following the adiabatic mapping calculations, structures along the approximate reaction path are examined carefully to ensure that the change of geometry is continuous along the path; subsequently, the saddle point is fully optimized by conjugated peak refinement (CPR)⁴³ to obtain more precise TS structure and energy. Finally, frequency calculations are carried out to confirm the nature of the stationary points and to compute the vibrational entropy and zero point energies to obtain approximate activation free energy; although using a harmonic approximation to estimate activation entropy is known to be of limited accuracy, previous studies of phosphate diester hydrolysis found that activation entropy does not differ much between different diesters²⁸. The vibrational frequencies are also used to estimate ¹⁸O kinetic isotope effects (KIEs) for MpNPP⁻ as an additional benchmark of the methodology (see Supporting Information).

To correct for intrinsic errors of SCC-DFTBPR energies, we explore corrections based on gas phase single-point energy calculations with both B3LYP/6-311++G** and MP2/6-311++G** at SCC-DFTBPR geometries; the B3LYP level was found to give very similar results to MP2 for simple phosphate hydrolysis reactions³⁷ (however, see below). As discussed in the literature,⁴⁴ such a simple correction may not always improve the energetics for semi-empirical methods given the errors in geometry; however, our previous tests^37,41 indicated that this correction scheme appears useful for SCC-DFTBPR since the method gives fairly reliable structures, even for transition states.

To further facilitate the analysis of sources of errors in both QM and QM/MM calculations, additional analysis for gas-phase/solution proton affinities (PA) for the leaving groups in MpNPP⁻ and its analogs. QM-only calculations are carried out by Gaussian03⁴⁵; PCM⁴⁶ and SM6⁴² models are employed to describe solvation effects. To test the accuracy of QM/MM coupling, solution PAs are also calculated using SCC-DFTB/MM based free energy perturbation calculations⁴⁷.

B. Enzyme Model Setup

For the hydrolytic reaction catalyzed by the AP superfamily members, a two-step mechanism is usually followed,⁴⁸ in which an oxygen nucleophile (e.g., Ser or Thr) first attacks the phosphorus/sulfur, then a water (hydroxide) replaces the leaving group in a step that is essentially the reverse of the first; for some family members of the superfamily, however, the mechanism can be more complex⁴⁹. In this work, to understand the catalytic mechanism of AP with phosphate diesters, we investigate the first step of the hydrolysis reaction of MpNPP⁻ in an E. coli AP variant in which Arg166 is mutated to Ser; this is expected to be the rate-limiting step given the experimental leaving group LFER analysis. Experimentally, this mutant was used to avoid inhibition by P_i,²² and it is believed that the mutation does not alter the reaction mechanism of AP since LFERs are similar in the mutant and WT.^20,50 Moreover, the chemical step is fully rate-limiting in this mutant. We also study a double mutant, R166S/E322Y AP, which was constructed in recent experimental studies to analyze the contribution(s) from the Mg²⁺ site in AP. For NPP, we study the wild type enzyme from Xanthomonas axonopodis pv. citri (Xac).

The enzyme models are constructed based on the X-ray structures for the E. coli AP mutant R166S with bound inorganic phosphate at 2.05 Å resolution (PDB code 3CMR³⁰) and Xac NPP with bound Adenosine Mono-Phosphate (AMP) at 2.00 Å resolution (PDB code 2GSU¹⁶). The ezyme model for R166S/E322Y AP is constructed based on the crystal structure of E322Y AP (PDB code 3DYC¹⁴) by mutating Arg166 to serine. In each case, starting from the PDB structure, the ligand is first “mutated” to the substrate of interest, MpNPP⁻ or MpNPPS⁻; two possible orientations of the substrates are considered for the AP active site, with the -OMe group oriented towards either the magnesium ion or Ser102 backbone amide (see additional discussions in Sect.III). Hydrogen atoms are added by the HBUILD module⁵¹ in CHARMM.⁵² All basic and acidic amino acids are kept in their physiological protonation states except for Ser102 and Thr90 in AP and NPP, respectively, which are assumed to be the neucleophiles and deprotonated in the reactive complex. Water molecules are added following the standard protocol of superimposing the system with a water droplet of 27 Å radius centered at Zn1²⁺ (see Fig.1 for atomic labels) and removing water molecules within 2.8 Å from any atoms resolved in the crystal structure.⁵³ Protein atoms in the MM region are described by the all-atom CHARMM force field for proteins⁵⁴ and water molecules are described with the TIP3P model.⁵⁵ The QM region includes groups most relevant to the reaction: the two zinc ions and their 6 ligands (Asp51, Asp369, His370, Asp327, His412, His331), Ser102 and MpNPP⁻ for R166S AP; for NPP, this includes two zinc ions and their 6 ligands (Asp54, Asp257, His258, Asp210, His363, His214), Thr90 and MpNPP⁻. Only side chains of protein residues are included in the QM region and link atoms are added between C_α and C_β atoms. A larger QM region also has been tested for R166S AP which further incorporates the entire magnesium site, including Mg²⁺, sidechains of Thr155, Glu322 and three ligand water molecules. Comparison of optimized structures using different QM regions indicates fairly similar optimized structures (see Supporting Information for details), thus the smaller QM region is used for the majority of the calculations. The treatment of the QM/MM frontier follows the DIV scheme in CHARMM; previous benchmark calculations have shown that this scheme generally gives reliable results for structure and energetics in QM/MM calculations provided that the MM charge is small near the QM/MM boundary⁵⁶ Since the Mg²⁺ near the QM region in AP is treated as a point charge (otherwise the QM region will become substantially larger), to avoid over-polarization of nearby QM groups, a NOE potential is added to the C-O bonds in Asp51, which is coordinated to both Mg²⁺ and Zn²⁺. The NOE potential takes the form:

in which R_min and R_max set the interval between which the restraining potential is zero; they are taken to be 0 and 1.28 Å, respectively. K_max is set to be 10⁴κcal/(mol · Ų).

Due to the fairly large size of the QM region (more than 80 atoms) and extensive sampling required for the open active site of AP and NPP, the SCC-DFTBPR method³⁷ is used for PMF calculations. Extensive benchmark calculations and applications indicate that it is comparable to the best semi-empirical method available in the literature for phosphate chemistry.^27,57

The generalized solvent boundary potential (GSBP)^58,59 is used to treat long range electrostatic interactions in geometry optimizations and MD simulations. The system is partitioned into a 27-Å spherical inner region centered at the Zn1 atom, with the rest in the outer region. Newtonian equations-of-motion are solved for the MD region (within 23 Å), and Langevin equations-of-motion are solved for the buffer regions (23-27 Å) with a temperature bath of 300 K; protein atoms in the buffer region are harmonically constrained with force constants determined from the crystallographic B-factors.⁶⁰ All bonds involving hydrogen are constrained using the SHAKE algorithm,⁶¹ and the time step is set to 1 fs. All water molecules in the inner region are subject to a weak GEO type of restraining potential to keep them inside the inner sphere with the MMFP module of CHARMM. The static field due to outer-region atoms, ${\varphi }_s^{io}$, is evaluated with the linear Poisson-Boltzmann (PB) equation using a focusing scheme with a coarse cubic grid of 1.2 Å spacing, and a fine grid of 0.4 Å spacing. The reaction field matrix M is evaluated using 400 spherical harmonics. In the PB calculations, the protein dielectric constant of ∊_p = 1, the water dielectric constant of ∊_ω = 80, and 0.0 M salt concentration are used; the value of ∊_p is not expected to make a large difference in this particular case because the active site is already very solvent accessible and the inner/outer boundary is far from the site of interest. The optimized radii of Nina et al.^62,63 based on experimental solvation free energies of small molecules as well as the calculated interaction energy with explicit water molecules are adopted to define the solvent-solute dielectric boundary. To be consistent with the GSBP protocol, the extended electrostatic model⁶⁴ is used to treat the electrostatic interactions among inner region atoms in which interactions beyond 12 Å are treated with multipolar expansions, including the dipolar and quadrupolar terms.

C. Benchmark enzyme calculations based on minimizations and reaction path calculations

To further test the applicability of SCC-DFTBPR/MM to AP and NPP, geometry optimization for the reactant (Michaelis) complex is compared to results from B3LYP^65–67/MM calculations. The basis set used in the B3LYP/MM calculations is 6-31G*⁶⁸, and the calculations are carried out with the QChem⁶⁹ program interfaced with CHARMM(c36a2 version).⁷⁰ Due to the rather large size of the QM region and the high cost of ab initio QM/MM calculations, atoms beyond 7 Å away from Zn1 are fixed to their crystal positions in these minimizations (note that these are not fixed in the potential of mean force simulations, see below. Also, test calculations at the SCC-DFTBPR level show that fixing atoms beyond 7 Å from Zn1 in minimizations do not lead to much difference as compared to a fully flexible inner-region calculation within the GSBP framework). The convergence criteria for geometry optimization are that the root-mean-square (RMS) force on mobile atoms is smaller than 0.30κcal/(mol · Å) and the maximum force smaller than 0.45κcal/(mol · Å).

The Minimum Energy Path (MEP) calculations are carried out by one-dimensional adiabatic mapping at both SCC-DFTBPR and B3LYP/6-31G* levels; the reaction coordinate is the antisymmetric stretch involving the breaking and forming P-O bonds (PO^lg-PO^nu), and the step size for the adiabatic mapping is 0.2 Å. At the SCC-DFTBPR level, the transition state is further refined using CPR.

D. 1D and 2D Potential of mean force (PMF) simulations

To study the free energy profile of enzyme reactions, PMF simulations have been carried out for R166S AP, R166S/E322Y AP and NPP with MpNPP⁻ and MpNPPS⁻ as the substrates. After the initial minimizations starting from the relevant crystal structure, the enzyme system is slowly heated to 300 K and equilibrated for 100 ps. The reaction coordinate is defined as PO^lg-PO^nu. The umbrella sampling approach⁷¹ is used to constrain the system along the reaction coordinate by using a force constant of 150 kcal/mol·Å⁻². In total, more than 51 windows are used for each PMF and 100 ps simulations are performed for each window. The first 50 ps trajectories are discarded and only the last 50 ps are used for data analysis. Convergence of the PMF is monitored by examining the overlap of reaction coordinate distributions sampled in different windows and by evaluating the effect of leaving out segments of trajectories. The probability distributions are combined together by the weighted histogram analysis method (WHAM)⁷² to obtain the PMF along the reaction coordinate. The averaged key structural properties for each window are calculated and summarized in Table IV.

In a separate set of PMF calculations, the Zn²⁺-Zn²⁺ distance is constrained to be 3.6, 4.1 and 4.6 Å, respectively, by a strong constraint with a force constant of 2,000 kcal/mol·Å⁻², to investigate the impact of this fundamental variable of the bimetallic site on catalysis in AP and NPP. For reference, the Zn²⁺-Zn²⁺ distance found in the various crystal structures for AP and NPP is close to be 4.1 Å.

To verify that the 1D PMFs capture the nature of the phosphoryl transfer transition state, we also carry out 2D PMF calculations for the α orientation of MpNPP⁻. The reaction coordinates are defined as the P-O^lg and P-O^nu distances, and the range of each distance is similar to that in the 1D PMF calculations. In total, 272 windows are used, and a force constant of 200 kcal/mol·Å⁻² is used for all windows; for each window, 100 ps simulations are carried out and only the last 50 ps are used for the subsequent WHAM analysis.

A. MpNPP⁻ hydrolysis in solution

The hydrolysis of MpNPP⁻ has been studied extensively by experiments and computations. Experimental studies²² determined the activation free energy of 25.7 kcal/mol at 42°C with hydroxide as the nucleophile, and the mechanism is established as concerted with a synchronous TS based on LFER analysis. Several computational work also studied the same reaction by employing various levels of theory. By using B3LYP/6-31+G* and the PCM model, Rosta and coworkers obtained a fairly loose TS with P-O^lg and P-O^nu as 1.86 and 2.49 Å, respectively.²⁸ By using B3LYP/6-311+G** with the COSMO continuum model on PCM-minimized geometries and a careful treatment of solute configurational entropy, they obtained an activation free energy barrier of 24.4 kcal/mol. In the more recent QM/MM simulations using explicit solvent (TIP3P) and AM1(d)-PhoT²⁷ as QM, López-Canut et al. obtained a free energy barrier of 20.5 kcal/mol; the transition state was featured with the P-O^lg and P-O^nu distances of 1.81 and 2.23 Å, respectively, somewhat more compact compared with the PCM result.

With our SCC-DFTB(PR)/PB method and charge dependent atomic radii⁴¹, the adiabatic map for MpNPP⁻ hydrolysis (Fig. 3a) is qualitatively consistent with previous studies and indicates a synchronous TS with an energy barrier of around 30 kcal/mol. After adding higher level (B3LYP or MP2) single point energy corrections, the general landscape of the adiabatic map does not change (Fig. 3b); the synchronous TS is still preferred with the barrier lowered to around 25 kcal/mol when MP2 corrections are used. After CPR refinement, the fully optimized TS (Fig. 3c) has a P-O^lg of 2.23 Å and P-O^nu of 2.43 Å; i.e., the P-O^lg distance is longer compared with previous theoretical results, while P-O^nu is similar. We note, however, the PES is rather flat near the transition state. The free energy barrier by including ZPE and solute configurational entropy is 29.3 kcal/mol at the SCC-DFTBPR level, which is decreased to 24.4 kcal/mol by including gas phase single point energy correction at the MP2/6-311++G** level (Table I); the latter value compares favorably with experimental value.

In addition, we study the hydrolysis of two other related diesters (Fig. 2), methyl 3-nitrophenyl phosphate (MmNPP⁻) and methyl phenyl phosphate(MPP⁻) with the approach. Experimentally, the trend is that the hydrolysis barrier increases as the pK_a of the leaving group increases (Table I)²². This trend has been reproduced by a previous theoretical study²⁸ in which the nature of transition states is found to be synchronous and becomes looser as the pK_a of the leaving group decreases; P-O^lg ranged from 1.84 to 1.86 Å and P-O^nu ranged from 2.33 to 2.49 Å. With our SCC-DFTBPR/PB method, this trend is also qualitatively reproduced, regardless of whether the gas-phase correction at higher level (B3LYP or MP2) is included. At a quantitative level, however, the SCC-DFTBPR/PB barriers are too high and the effects of the substitution are overestimated (see Table I); for example, the barrier difference between MpNPP⁻ and MPP⁻ is only 2.9 kcal/mol according to experiment, but 10.5 kcal/mol at the SCC-DFTBPR/PB level. The discrepancy remains fairly large even with B3LYP corrections, while including the MP2 gas-phase corrections significantly improves the agreement with experimental value. The nature of the TS also becomes somewhat tighter (especially P-O^lg, by ∼ 0.2 Å) when MP2 correction is included (see Supporting Information).

To better understand the quantitative differences between SCC-DFTBPR, B3LYP and MP2 results, we examine the relative PAs of the leaving groups in the three phosphate diesters in both gas-phase and solution; gas-phase PAs reflect the intrinsic accuracy of the QM method, while solution PA calculations also examine the accuracy of either the implicit solvent model or QM/MM interactions in explicit solvent simulations⁷³. As shown in Table II, all DFT methods, which include both SCC-DFTB(PR) and B3LYP, have errors much larger than “chemical accuracy” (1 kcal/mol) for the relative gas-phase PAs, especially concerning the effect of introducing the nitro group; by comparison, MP2 does a much better job. The errors in the relative solution PAs follow the same trend as the relative gas-phase PAs, suggesting that errors in the gas-phase PAs are the major source of error; this is confirmed by the observation that computed relative solvation free energies for the leaving groups considered here are in good agreement with experimental values using both SCC-DFTB(PR) and B3LYP based implicit solvent models, with the exception of B3LYP and UAKS radii (see Supporting Information).

As additional benchmark for the nature of the TS, ¹⁸O KIE calculations for MpNPP⁻ in solution are carried out with the SCC-DFTBPR transition state and harmonic vibrational frequencies. As shown in Table S2, the trends are in qualitative agreement with experimental results^74,75 and previous AM1(d)-PhoT and B3LYP results.⁷⁶ On the other hand, we note that SCC-DFTBPR overestimates the magnitude of the KIEs, especially for the effects associated with the leaving group oxygen and the non-bridging oxygen. This is consistent with the trend discussed above that SCC-DFTBPR predicts a solution TS that is looser as compared to previous theoretical calculations, with most notably a weaker and longer P-O^lg bond.

In short, the benchmark calculations for the phosphate diesters and their leaving groups suggest that SCC-DFTBPR can provide fairly reliable structural properties of these species and a semi-quantitative description of energetics and the nature of hydrolysis transition state, especially for relative trends associated with different substituents on the leaving group.

B. First step of MpNPP⁻ hydrolysis in R166S AP

Based on the crystal structure of the AP R166S mutant complexed with inorganic phosphate, the phosphate ligand is “mutated” to MpNPP⁻ by adding necessary functional groups to phosphate oxygen. The leaving 4-nitrophenyl group is added to O1 due to the geometrical requirement of the in-line attack from Ser102 (see Fig. 1a). The methyl group can be added to O3 or O4, which correspond to two different substrate orientations (denoted as α and β orientations, respectively, following the notation of Ref.¹⁴). Recent experimental studies¹⁴ using a double mutant AP with the Mg²⁺ site removed and phosphorothioate diesters suggested that the α orientation is preferred over the β orientation. As discussed below (see Sect.III C), however, the interpretation of those elegant experiments may not be as clearcut as presented. Moreover, even if the α orientation is indeed dominant, it is not clear if the discrimination comes from binding or the chemical step. Therefore, it is informative to study both orientations.

The comparison of optimized structures by B3LYP/MM and SCC-DFTBPR/MM shows good agreement between the two levels (Fig.4a). The O_Ser102-P distance increases from 3.1 Å in the crystal structure to 3.4 (3.3) /3.8 (3.9) Å in B3LYP/MM (SCC-DFTBPR/MM) optimized structure for the α/β orientation, leading to a stable reactant complex. The O2 of the substrate coordinates to one of the zinc ions and O1 with the phenyl group is solvated by water molecules. Interestingly, the slight shift of the substrate position also increases the distance between a Mg²⁺-bound water (Wat1) and substrate O3 such that a hydrogen bond is formed between Wat1 and the Zn²⁺-activated Ser102, instead of with an inorganic phosphate oxygen as in the crystal structure (this holds also in MD simulations at the SCC-DFTBPR/MM level; also see discussion below for interactions in the TS). O4 and the nearby Ser102 backbone amide forms the only direct hydrogen bond between the substrate and the enzyme, which is shorter (1.9 Å vs. 3.3 Å) with the substrate in the α orientation than in the β orientation. If this is the only major difference between those two orientations, the binding affinity of MpNPP⁻ to the enzyme is likely stronger with the α orientation than with the β orientation, although a more quantitative estimate remains to be carried out with free energy simulations, which we defer to a future study. Many other hydrogen-bonding distances (e.g., between Wat1 and Ser102, MpNPP⁻ and Ser102 backbone amide) and distances involving the zinc ions (e.g., between Ser102/MpNPP⁻ and the zinc ions) are similar at the two levels of theory (see Fig.4). The Zn²⁺-Zn²⁺ distance is generally shorter at the SCC-DFTBPR/MM level (by ∼0.1-0.2 Å) while the distances between Zn and its ligand oxygen are generally longer. Overall, however, the agreement between optimized structures at the two levels of theory is excellent, supporting the use of SCC-DFTBPR/MM.

In addition to the structural similarity in optimized reactants, the MEP results (for β orientation) from adiabatic mapping also show good agreement (16.9 vs. 15.7 kcal/mol) between SCC-DFTBPR/MM and B3LYP/6-31G*/MM calculations, which further supports the use of SCC-DFTBPR/MM; the adiabatic mapping and CPR calculations at the SCC-DFTBPR/MM level give similar transition states (see Supporting Information), although the barrier height is slightly lower with fully (CPR) optimized saddle point (see Table III). As shown in Fig.4c, the main differences between SCC-DFTBPR/MM and B3LYP/MM transition states include a shorter Zn²⁺-Zn²⁺ distance at the former level, a shorter P-O_Ser102 bond length and hydrogen-bond distances for the interaction between the substrate and nearby groups (e.g., Ser102 backbone amide and Mg²⁺-bound Wat1). We note that, in those MEP calculations with SCC-DFTBPR, the Zn²⁺-Zn²⁺ distance appears somewhat shorter in the TS for the β orientation than the α orientation (see Supporting Information), although the difference is much smaller in the PMF calculations (see below), again highlighting the importance of sampling protein fluctuations. Nevertheless, there is room to further improve the SCC-DFTBPR method, which may require including complete third-order terms in the SCC-DFTB expansion⁷⁷ and a more systematic refitting of the P-O repulsive potential⁷⁸.

For both substrate orientations, the PMF peaks at the reaction coordinate (PO^lg-PO^nu) slightly less than 0 Å and then drops in the product region, corresponding to an exothermic reaction (Fig. 5a, 7a). The free energy barriers are 23.4 and 19.6 kcal/mol for the α and β orientations, respectively (see Table III). It is worth mentioning that these barriers correspond to the free energy difference between the TS and the Michaelis complex, i.e., κ_cat, while experimentally reported values are κ_cat/K_M, which prevents a direct comparison between calculation and experiment. Nevertheless, the measured κ_cat/K_M for R166S AP with MpNPP⁻ corresponds to a free energy barrier of 18.0 kcal/mol, which gives the lower bound for the free energy barrier for the chemical step. Therefore, our calculated barriers for the chemical step are qualitatively consistent with the measured κ_cat/K_M value. The calculations also suggest that the α orientation has a higher barrier for the chemical step. Therefore, for the α orientation to be at least as competitive as the β orientation in terms of κ_cat/K_M, the corresponding binding free energy should be at least 3.8 kcal/mol stronger, which is qualitatively consistent with the above observation of a stronger hydrogen bonding interaction between the enzyme and MpNPP⁻ in the α orientation. Although further binding free energy calculations need to be carried out, the results suggest that the model¹⁴ in which the α orientation is the only productive binding mode seems oversimplified (see Sect.III C below for additional discussions).

The key structural properties of the active site are averaged over the trajectory of each window and plotted as functions of the reaction coordinate (Fig. 5b, 7b). The changes of P-O^lg and P-O^nu clearly show that the concerted pathway with a synchronous TS is operative for both substrate orientations and similar to the reaction in aqueous solution, supporting that AP does not significantly alter the nature of hydrolysis TS for phosphate diesters; this is in qualitative agreement with LFER data for R166S AP and a series of substituted methyl phenyl phosphate diesters²². The nature of the calculated TS is not sensitive to the definition of the 1D reaction coordinate, since the 1D (Fig.5) and 2D (Fig.6) PMF calculations show very consistent results in terms of both structures and energetics. In terms of the tightness coordinate (TC=PO^lg+PO^nu), it decreases from 4.66 Å in solution to 3.89 (α) and 3.94 Å (β) in R166S AP (Table IV); this decrease is likely due to the bi-metallic zinc motif in AP, since both MpNPP⁻ and Ser102 are coordinated with Zn²⁺ ions, which are separated by ∼ 4 Å throughout the reaction. We note that the degree of tightening is likely not as large as the values for TC imply since our calculations appear to overestimate the value of TC for solution reactions as compared to previous calculations^25,28.

We note that recent QM/MM simulation of WT AP with both mono- and di-esters^24,26 showed rather large structural changes compared to the crystal structure; for example, the Zn²⁺-Zn²⁺ distance increases up to 7 Å, while no such major distortion is observed here (also see discussions below in Sect.III E). In the TS for the β orientation (Fig.7d), Wat1 of the Mg²⁺ breaks the hydrogen bond with Ser102 and forms a new one with O4 of MpNPP⁻, which is presumably due to the larger reduction of charge on the Ser102-O than the substrate O4; the average Mulliken charges for Ser102 O and O4 are -0.79 and -0.92, respectively, in the reactant, -0.56 and -0.89, respectively, in the TS. This change of hydrogen-bonding interactions likely helps lower the reaction barrier compared with the α orientation in which Wat1 interacts loosely with both Ser102-O and O3 in the TS (see Fig.5d); the average Mulliken charges for Ser102 O and O3 are -0.79 and -0.32, respectively, in the reactant, -0.50 and -0.40, respectively, in the TS.

Another interesting observation is that the leaving group does not form a direct interaction with any zinc ion in either the reactant state or TS; rather, it is “solvated” by water molecules accessible to the fairly open active site. This binding mode (especially in the TS) is in contrast to that observed for vanadate, a widely used transition state analog for phosphoryl transfers, in the crystal structures for AP-vanadate and NPP-vanadate complexes^16,79; these structures suggest a binding mode in which one non-bridging oxygen and the leaving oxygen interact directly with one of the zinc ions. Benchmark calculations suggest that our QM/MM protocol is able to reproduce the binding mode of vanadate in the active site of both AP and NPP, and that SCC-DFTBPR describes the interaction between the di-metallic zinc motif and phosphate diesters in good agreement with B3LYP (see Supporting Information). Moreover, MD simulations starting from the vanadate binding mode with the reaction coordinate constrained to be zero converge to the same binding mode shown in Figs.5-7. Collectively, these results suggest that the binding mode observed in the current work is unlikely an artifact of the computational methodology and indeed energetically favorable for systems studied here. We note that the TS for diesters in AP (and NPP, see below) is rather tight in nature, thus the leaving group oxygen doesn't bear any significant formal charge. Therefore, the binding mode captured in the vanadate structures better reflects the situation for monoesters, which feature a much looser TS in which the leaving group is substantially more charged. Our preliminary calculations for monoesters in AP indeed find tighter interactions between the zinc ion and the leaving group (Hou and Cui, work in progress).

C. Additional analysis of substrate orientation: activity in the double mutant (R166S/E322Y) and thio effects in R166S AP

To further clarify the issue of substrate orientation in AP, we carry out simulation studies to analyze two sets of experiments that were designed to answer the same question. In the first set of experiments, Zalatan and coworkers constructed mutants with the Mg²⁺ site removed (E322Y, E322A, R166S/E322Y) and measured the catalytic activities for phosphate monoester, phosphate diester and sulphate monoester in these mutants¹⁴. Based on the observed large detrimental effects of Mg²⁺ removal on phosphate monoester and sulphate monoester hydrolysis but negligible effect on phosphate diester hydrolysis, they concluded that the α orientation is preferred over the β orientation (also see the discussion on thio effects below). However, these observations alone do not rule out the possibility that the two orientations are in fact similar in activity in the WT (and R166S) enzyme, which is what we have observed in this study (see Table III). In a mutant where one pathway is significantly perturbed, the other can still provide an alternative route, which may explain the only 2 fold decrease in κ_cat/K_M for R166S/E322Y AP as compared to R166S AP (both with MpNPP⁻ as the substrate). To support this, we explicitly carry out calculations for the double mutants.

As mentioned in Computational Methods, the R166S/E322Y double mutant simulations are prepared based on the crystal structure of E322Y AP¹⁴, in which the hydroxyl group of Tyr322 occupies the region corresponding to the Mg²⁺ site in WT AP and forms a hydrogen bond with Asp51; the bimetallic zinc site is largely unaffected. Based on the comparison of results for the α and β orientations in R166S AP (Figs.5-7), we expect that mutating away the Mg²⁺ site, which turns off interactions between the Mg²⁺-bound water and substrate oxygen, will result in a large detrimental effect on the β orientation but a much smaller effect on the α orientation. This is exactly what we observe from the double mutant calculations; as shown in Table III, the reaction barrier is slightly decreased to 22.6 kcal/mol for the α orientation but is increased to be over 30 kcal/mol for the β orientation. These results directly support the important role of the Mg²⁺ site in reducing the barrier for the β orientation, and by inference, for the hydrolysis of phosphate monoesters, as observed experimentally¹⁴. Analysis of structures from PMF calculations also indicates that a water molecule penetrates into the double mutant active site for the β orientation to further stabilize the nucleophile in the reactant state, thus also contributing to the significant increase of the barrier; no such water penetration is observed for the α orientation since the hydrophobic -OMe group in the substrate helps block additional water from the active site. The double mutant calculations also explicitly support that the α orientation is not affected much by the Mg²⁺ site, in agreement with the experimental observation¹⁴ that the activity of phosphate diesters remains largely unperturbed in the Mg²⁺-site mutants. In short, our model that both α and β orientations are productive binding modes in R166S AP while only α is reactive in the Mg²⁺-site mutants is consistent with available experimental data.

In the second set of experimental studies, Zalatan and coworkers analyzed the reactivities of phosphorothioate diesters in several variants of AP¹⁴. Take MpNPPS⁻ as an example (Fig.8), the key observation was that R166S AP reacts with the R_p enantiomer at least 10² times faster than with the S_p enantiomer, suggesting that the binding mode with the R′ (Me) group toward the Mg²⁺ site and the sulfur toward the Ser102 backbone amide is dominant; a relevant piece of information here is that previous experiments in the same group established that the sulfur is not placed between the Zn²⁺ ions¹⁵. We note, however, these discussions rely on the assumption that phosphorothioate diesters behave, in terms of hydrogen bonding with active site groups, similarly to phosphate diesters, which may not be as clearcut as commonly believed. For example, NBO charge analysis in both gas phase and solution (see Fig. 8) shows that in MpNPPS⁻ the oxygen bonded with the methyl group is in fact more negatively charged than the sulfur; i.e., it is not obvious that the R_p enantiomer of MpNPPS⁻ reflects the a-orientation of MpNPP⁻. Therefore, although the experimental observation that removal of the Mg²⁺ site (R166S/E322Y AP vs. R166S) has a similar impact on the hydrolysis of MpNPPS⁻ and MpNPP⁻ suggests that interactions from the Mg²⁺ site are similar in the phosphorothioate and phosphate esters, it is not clear if the thio effects can unambiguously infer the binding mode of phosphate diester.

To help better understand the thio effects, we calculate the PMF for MpNPPS⁻ hydrolysis in both R166S AP and R166S/E322Y AP. The results support that in R166S AP, the R_p enantiomer is indeed favored over S_p enantiomer by 7.9 kcal/mol for the chemical step (see Table III), in qualitative agreement with experimental findings. The experimental thio effect (the ratio of rate constant for phosphate ester substrate over that for the phosphorothioate analog) corresponds to a free energy barrier difference of 3.6 kcal/mol, while our calculated value is 2.1-5.9 kcal/mol, depending on whether the α or β orientation for MpNPP⁻ is used as reference; note again, however, the experimental value is based on κ_cat/K_M while our values are based on the chemical step only. One point worth mentioning is that due to the weaker substrate binding following thio substitution, the Solvent Accessible Surface Area (SASA) of the sulfur in MpNPPS⁻ is much larger than that for the corresponding oxygen in MpNPP⁻, especially for the S_p enantiomer (see Table 2 in Supporting Information); this does not occur at the transition state. The larger SASA provides extra solvent stabilization for the reactant state and probably accounts partially for the much larger thio effects calculated for the S_p enantiomer (13.8 kcal/mol higher in barrier relative to the β orientation of MpNPP⁻). The PMF profiles for MpNPPS⁻ also peak at where the value of the reaction coordinate is slightly less than 0 while the position of the reactant state is decreased from ∼2 Å in MpNPP⁻ to ∼-2.5 to -3 Å, reflecting the larger substrate size and weaker binding interactions with the active site than MpNPP⁻ (see Supporting Information); the transition state is still synchronous in nature and slightly tighter than MpNPP⁻.

The effect of removing the Mg²⁺ site (in R166S/E322Y AP) on the hydrolysis of MpNPPS⁻ is expected to be small for both enantiomers since the interactions between Mg²⁺-water and either the -OMe or -S⁻ are fairly weak. Indeed, for the R_p enantiomer, the barrier actually decreases in the R166S/E322Y mutant relative to the R166S AP to 19.7 kcal/mol; this value is slightly lower than the experimental κ_cat/K_M value of 21.0 kcal/mol. For the S_p enantiomer, surprisingly, the barrier increases to be over 40 kcal/mol. A closer examination of the simulation snapshots shows that for the reactant state with the S_p enantiomer a solvent water penetrates into the active site and forms a hydrogen bond with the nucleophile (Ser102 oxygen), as illustrated by the comparison of integrated water distribution near the nucleophilic oxygen in Ser102 in the reactant state (see Supporting Information). As mentioned above for the β orientation of MpNPP⁻ in R166S/E322Y AP, which also features a significantly increased barrier (Table III), a similar water penetration is observed as well. The water penetration to Ser102 only happens for these two cases and is reproducible in simulations in which the penetrated water is first deleted and then the system further equilibrated. We note that in R166S AP the active site is already rather open to solvent molecules, thus additional water penetration into the active site is not unexpected when the Mg²⁺ site is removed; the β orientation and the S_p enantiomer are particularly susceptible to water penetration since they lack the bulky methyl group near Ser102. Nevertheless, water penetration in AP mutants remains an interesting issue that deserves in-depth analysis from future experimental and computational studies.

Considering the results for MpNPP⁻ and those for the R_p enantiomer for MpNPPS⁻ in the two AP variants, our calculations qualitatively reproduced key experimental observations concerning the effects of Mg²⁺ site removal and thio substitution, further supporting the argument in the last subsection that experimental data so far can not be used to unambiguously determine the orientation of diester substrates in the AP active site. In the broader context, as we mentioned above, the charge distributions for MpNPP⁻ and MpNPPS⁻ bear some nontrivial differences (Fig.8); in addition, the possibilities of water penetrating into the active site for certain orientation of the (thio-substituted) substrate and that different substrate orientations are dominant in different variants of AP (R116S vs. R116/E332Y) further complicate interpretation of the observed thio effects.

D. First step of MpNPP⁻ hydrolysis reaction in NPP

The hydrophobic groove in NPP has been suggested to contribute at least 10⁴-fold to the catalysis of phosphate diester reactions¹⁶ by favorable interactions with the extra R′ group in diesters (Fig.1b). Therefore, only one orientation is studied here for MpNPP⁻ in which the methyl group points toward the hydrophobic pocket. Similar to the comparisons made above for AP, SCC-DFTBPR/MM minimizations for MpNPP⁻-NPP also give similar reactant complex structure to B3LYP/MM calculations (Fig.4b). The O_Thr90-P distance increases from 3.2 Å in the crystal, which contains AMP as the inhibitor, to 3.6 (3.6) Å at the SCC-DFTBPR/MM (B3LYP/MM) level. The substrate O2 coordinates with Zn1 while O4 forms hydrogen bonds with Asn111 and the backbone amide of Thr90. The optimized Zn²⁺-Zn²⁺ distance is 4.47 (4.40) at the B3LYP/MM (SCC-DFTBPR/MM) level. The two hydrogen bonds formed between O4-Asn111 and O4-Thr90-backbone-amide are also in decent agreement at different levels of theory (Fig.4b).

The PMF calculation shows that the free energy profile corresponds to an exothermic process with the barrier located at PO^lg-PO^nu ∼-0.20 Å and a barrier height of 20.2 kcal/mol. The measured κ_cat/K_M corresponds to 14.3 kcal/mol at 300 K¹⁶ and sets the lower limit for the chemical step barrier. Compared with AP, the calculated barrier for NPP is close to the β orientation but lower than that in the α orientation. Since MpNPP⁻ is the cognate substrate of NPP and therefore expected to bind tighter to NPP than to AP, the calculated barrier implies a higher κ_cat/K_M value for NPP than for AP, which is consistent with experimental observations. In other words, the calculations explicitly support that although diesters are cognate substrates for NPP and promiscuous substrate for AP, the chemical step in NPP is not much accelerated over (R166S) AP.

The changes of P-O^lg and P-O^nu (Fig. 9b) show that the concerted pathway is also operative for NPP, with a TS similar to that in aqueous solution and AP. For example, the tightness coordinate is 3.86 Å (Table IV), as compared to the values of 3.89(α)/3.94(β) and 4.66 Å, respectively, in R166S AP and solution, respectively. Considering that the tightness coordinate for solution TS seems overestimated by our method compared to previous work^25,28, our calculations support the idea motivated by LFER data²² that, instead of altering transition state structure, NPP and AP catalyze phosphoryl transfer reactions by recognizing and stabilizing transition states similar to those in aqueous solution.

E. Comparison to recent QM/MM simulations^25,26

As just stated, our calculations find that the transition states for diester hydrolysis in AP and NPP are similar and slightly tighter than that in solution. This is in direct contrast to the recent QM/MM studies^25,26 which found that the TS in AP/NPP is much looser in nature with the tightness coordinate of 5.66 (NPP)/5.00 (AP) vs. a value of 4.05 in solution. Several pieces of evidence suggest that those calculations are less reliable than our SCC-DFTBPR/MM calculations. First, as noted above for AP, their calculations led to large structural distortions in the bi-metallic zinc motif relative to the crystal structure, while no such distortions occur in our calculations; the same trends hold for NPP calculations, and the loose TS found in previous work^25,26 might be a result of the substantially elongated Zn²⁺-Zn²⁺ distance. Another example concerns the hydrogen-bond between Asn111 and MpNPP⁻ in NPP, which is observed in the crystal structure (with AMP and vanadate) and throughout our simulations; by contrast, this interaction was broken in the recent QM/MM simulations.²⁵ Second, our calculated barrier heights are consistently higher than the barriers that correspond to experimentally measured κ_cat/K_M values, while this is not the case for recent QM/MM calculations for both NPP²⁵ and AP²⁶. For example, for NPP with the same substrate, their best estimate of the barrier²⁵ was ∼11 kcal/mol, which is even lower than the barrier estimated based on experimental κ_cat/K_M, casting further doubt on the quantitative nature of their result.

What could be the origin for the differences between our and the recent QM/MM calculations^25,26? Both SCC-DFTBPR and AM1(d)-PhoT are approximate in nature and parameterized rather specifically for phosphate hydrolysis reactions. For solution reaction, the AM1(d)-PhoT approach gives encouraging results in terms of both energetics and KIEs⁷⁶, while SCC-DFTBPR, once corrected for gas-phase energies based on MP2, gives good barriers but too loose transition states as reflected by the computed KIE values. In the context of zinc containing enzymes, however, the SCC-DFTBPR approach has been tested with more benchmark calculations, including those discussed here. Indeed, although AM1/PM3 has been used successfully to describe a number of metalloenzymes that catalyze phosphoryl transfers^80,81, combining AM1 for zinc and AM1(d)-PhoT for phosphoryl transfers in zinc enzymes has not been carefully tested. Therefore, we suspect that the main cause for the large structural variations in the recent QM/MM simulations^25,26 is the use of AM1 for zinc. In this regard, although it is possible that the crystal structure with an inhibitor (e.g., inorganic phosphate) doesn't capture all structural features (e.g., variation in the zinc-zinc distance) of the transition state, the fact that the active site undergoes little change with a transition state analogue (vanadate, see discussions in Sect.III B and also Supporting Information) in both high-resolution crystallography and EXAFS studies³¹ suggests that large variations in the zinc-zinc distance seen in the recent QM/MM simulations^25,26 (including in the Michaelis complexes!) are unlikely realistic. Since the large increase in the zinc-zinc distance occurs in their calculations of several AP and NPP variants^25,26, the effect likely cancels out for relative trends; this explains why mutation effects were adequately captured in Ref.²⁶. Finally, we noted that the QM/MM boundary in Refs.^25,26 cuts across fairly polar covalent bonds yet the simple single link-atom scheme was used. As emphasized by several researchers^56,82–84, extra care needs to be exercised when the QM/MM boundary involves polar bonds. This technical detail likely also contributes to the uncertainty of the results in Refs.^25,26.

F. Why is the nature of TS for phosphate diesters in AP and NPP similar to that in solution?

As another way to evaluate the conflicting findings in the current and previous QM/MM studies, we further dissect whether our observation that the nature of TS for phosphate diesters in AP and NPP is similar to that in solution is consistent with other known experimental facts. We do so with the scheme outlined in Fig.11, which is qualitatively similar to that used by Herschlag and co-workers^18,22,23.

For aryl phosphate diester hydrolysis in solution, the measured reaction barrier, ΔG^‡(aq) is ∼26 kcal/mol (e.g., see Table I for MpNPP⁻); the corresponding barrier in the enzyme, ΔG^‡(E), is, for R166S AP, 18.0 kcal/mol (see κ_cat/K_M in Table III). Therefore, the enzyme binds to the TS, which is shown here to be of rather similar synchronous nature in solution and AP/NPP (Table IV), by about 8 kcal/mol $(\Delta \Delta G_{\mathit{\text{syn}}‡}^b)$.

For the enzyme to shift the nature of TS from synchronous to loose, the driving force needs to be large enough to overcome the binding energy for the synchronous TS $(\Delta \Delta G_{\mathit{\text{syn}}‡}^b)$ plus the energy gap between these two kinds of structures in solution $(\Delta \Delta G_{\mathit{\text{syn}}/\mathit{\text{loose}}}^{‡(aq)})$. The latter, although not measurable directly with experiments, can be estimated based on calculations; our calculations shown in Fig.3 give a value ∼ 8 kcal/mol (the “loose” structure is taken to have a tightness coordinate of ∼5.7Å as found for the TS in NPP in previous QM/MM calculations²⁵). In other words, the enzyme needs to bind to the loose TS by more than 16 kcal/mol $(\Delta \Delta G_{\mathit{\text{syn}}‡}^b+\Delta \Delta G_{\mathit{\text{syn}}/\mathit{\text{loose}}}^{‡(aq)}=8+8)$ to make the loose TS more favorable in the enzyme than a synchronous one.

Considering what we know about the activity of AP toward phosphate monoesters, however, we argue that such a strong binding is unlikely for phosphate diesters. For a phosphate monoester related to the diesters studied here, such as pNPP²⁻ (p-nitrophenyl phosphate), the solution barrier is about 32 kcal/mol¹⁸ and the barrier in R166S AP is 10.6 kcal/mol⁸⁵. Since LFER data indicate that the nature of TS is loose in both AP and solution^20,85, these results suggest that R166S AP binds a loose TS for monoesters by ∼21 kcal/mol. Since diesters feature less charge and are promiscuous substrates of (R166S) AP, we expect that the binding energy of a loose TS for diesters to R166S AP is substantially lower than 21 kcal/mol. Therefore, we don't expect that R166S AP is able to shift the TS for diesters to be much looser in nature, in agreement with our QM/MM calculations.

For NPP, the diesters are cognate substrates, thus it is conceivable that the active site has been evolved to optimize the catalysis of diester hydrolysis, which might involve modifying the nature of TS relative to solution. However, we note that at least for the diester substrate studied here, the calculated chemical step barrier in NPP (20.2 kcal/mol) is not significantly lower than that in solution (25.7 kcal/mol). Therefore, it is not unreasonable that the nature of TS in NPP is not significantly changed relative to solution and a significant component of rate enhancement over solution (κ_cat/K_M relative to κ_ω) is due to substrate binding.

G. The effects of Zn²⁺-Zn²⁺ distance on reaction energetics

Since the bimetallic zinc site is a prevalent catalytic motif^86–90, it is of interest to establish what features (e.g., structural vs. electrostatic) are important to the catalytic proficiency. In this work, motivated in part by the fact that the Zn²⁺-Zn²⁺ distances behave rather differently in the current and previous QM/MM calculations^24,25 of AP and NPP, we examine the effect of this fundamental structural feature on the catalysis in AP and NPP.

Specifically, we design two sets of simulations for R166S AP and NPP with MpNPP⁻ (α orientation only) to study the effects of Zn²⁺-Zn²⁺ distance fluctuation and variation. In the first set, we constrain the Zn²⁺-Zn²⁺ distance close to the value in crystal structures of AP and NPP (4.1 Å). Due to this constraint, the Zn²⁺-Zn²⁺ distance fluctuations in PMF simulations are significantly damped (root-mean-square-fluctuation, rmsf, of ∼0.01 Å) as compared to unconstrained simulations (rmsf∼0.2 Å). As shown in Fig.10a-b, the PMFs for AP and NPP do not exhibit much change, especially in the barrier height, due to the constraint; the nature of the TS (see Supporting Information) also does not change. These observations suggest that fluctuation of the Zn²⁺-Zn²⁺ distance during the reaction is not critical to the reaction barrier or the nature of the TS. The significant overlaps of the PMFs also indirectly support the reproducibility and convergence of our PMF simulations.

In the second set of simulations, the Zn²⁺-Zn²⁺ distances are constrained at 3.6 and 4.6 Å, respectively, which represent the two extreme values observed in unconstrained simulations. The changes of PMF are similar in AP and NPP (see Fig.10c-d): the reactant state position is shifted to a less negative value and the barrier height is reduced as the Zn²⁺-Zn²⁺ distance is decreased from 4.6 to 3.6 Å. In terms of the nature of TS, there is no qualitative change as reflected by the essentially invariant peak position of the PMFs. At the quantitative level, there are variations as reflected by the average tightness coordinate (see Table IV) and other structural parameters (see Supporting Information); for example, the tightness coordinate (TC) increases, as expected, as the Zn²⁺-Zn²⁺ distance increases.

In short, these calculations explicitly demonstrate that the Zn²⁺-Zn²⁺ distance of the bimetallic zinc site plays an important role in tuning the catalysis. This is not unexpected since the distance between the zinc ions influences the electrostatic properties in the active site, which are crucial to the phosphoryl transfers⁹¹. Nevertheless, the results clearly underlines the importance of reproducing geometrical properties of a bimetallic site for a meaningful analysis of its catalytic properties.

H. Issues worthwhile investigating with future experiments

The key objective of this work is to characterize the nature of the hydrolysis TS of phosphate diesters in AP and NPP so as to evaluate different hypotheses^22,26 regarding the catalytic promiscuity in these enzymes. To evaluate the main findings of this work, we propose that the following experimental studies are worthwhile.

First, unlike the recent QM/MM simulations^25,26, our calculations suggest that the nature of diester hydrolysis TS is largely similar in AP, NPP and solution. Along this line, LFER and KIE studies for diester hydrolysis in NPP will be highly informative.

Second, our calculations find that the leaving group does not interact strongly with the zinc ion in either the reactant or TS for diester hydrolysis in AP/NPP. This is rather unexpected considering the crystal structures of vanadate-AP/NPP complexes^16,79 and therefore worth further investigations. Note that this result does not suggest that the second zinc ion does not play an essential role in catalysis since it coordinates with the bridging oxygen. The leaving group is stabilized by active site solvent molecules, suggesting that the dependence of the catalytic rate with respect to the substitution of the leaving group, including both chemical and isotope substitutions, should be close to that in solution.

We have extensively discussed in Sect.III C the issue of different binding orientations of the substrate and relation of current calculations with available experimental studies. As highlighted by the NBO charges shown in Fig.8, one complicating factor for the interpretation of thiol experiments is that phosphorothioate and phosphate esters have rather different charge distributions; in MpNPPS⁻ the oxygen bonded with the methyl group is in fact more negatively charged than the sulfur. Therefore, it will be valuable to explore the reactivity of thiol substrates with both well-defined stereochemistry and relative charge distributions closer to the phosphate esters, such as by replacing both O3 and O4 in MpNPP by sulfur.

Finally, as alluded to in Sect.III G, the metal-metal distance appears to be essential to both the barrier height and the nature of the transition state; the barrier is higher and the nature of TS looser with a longer metal-metal distance. This can be tested by substituting the zinc ions with other metals of different size (e.g., replacing Zn²⁺ by Co²⁺) and performing kinetic and LFER analyses. Along this line, analysis of TS nature with synthetic analogs of di-metallic zinc motifs^92,93 will be highly informative, since the structure of such artificial catalysts can be better controlled and monitored.