|Home | About | Journals | Submit | Contact Us | Français|
Combined QM(PM3)/MM molecular dynamics simulations together with QM(DFT)/MM optimizations for key configurations have been performed to elucidate the enzymatic catalysis mechanism on the detoxification of paraoxon by phosphotriesterase (PTE). In the simulations, the PM3 parameters for the phosphorous atom were re-optimized. The equilibrated configuration of the enzyme/substrate complex showed that paraoxon can strongly bind to the more solvent-exposed metal ion Znβ, but the free energy profile along the binding path demonstrated that the binding is thermodynamically unfavorable. This explains why the crystal structures of PTE with substrate analogues often exhibit long distances between the phosphoral oxygen and Znβ. The subsequent SN2 reaction plays the key role in the whole process, but controversies exists over the identity of the nucleophilic species which could be either a hydroxide ion terminally coordinated to Znα or the μ-hydroxo bridge between the α- and β-metals. Our simulations supported the latter and showed that the rate-limiting step is the distortion of the bound paraoxon in order to approach the bridging hydroxide. After this preparation step, the bridging hydroxide ion attacks the phosphorous center and replaces the diethyl phosphate with a low barrier. Thus, a plausible way to engineer PTE with enhanced catalytic activity is to stabilize the deformed paraoxon. Conformational analyses indicate that Trp131 is the closest residue to the phosphoryl oxygen, and mutations to Arg or Gln or even Lys which can shorten the hydrogen bond distance with the phosphoryl oxygen could potentially lead to a mutant with enhanced activity for the detoxification of organophosphates.
Pesticides have been extensively used to control the spread of unwanted insects or weeds and generally categorized depending on their specific functions. Among the various groups of pesticides, organophosphates are one of the most widely used types, accounting for about 36% of the total world market, with more than 100 compounds being commonly used. Their high effectiveness has led to worldwide use, which has subsequently caused serious environmental problems.1 Although chemical structures differ within or between pesticide categories, they share a common undesirable feature, namely that they are often inherently toxic to non-target organisms, including human beings. Among the most popular organophosphorus pesticides, paraoxon and parathion are chemical nerve agents which inhibit mammalian acetylcholinesterase (AChE) through the irreversible phosphorylation of an active site serine residue. The effects associated with these chemicals may be either immediate or delayed.2 Organophosphorus esters such as tabun, sarin and O-ethyl S-(2-diisopropyl aminoethyl) methylphosphonothioate (VX in short) which are forerunners of the organophosphorus pesticides, however, are also being used as warfare nerve gases.3,4 Since toxic pesticide waste and chemical stockpiles pose a serious potential threat to both the environment and human health,5 it is of great interest and significant importance to develop effective and economical methods for the detoxification and removal of these organophosphates from nature and the battlefield.
Biodegradation provides a safe and efficient way to clean the residual organophosphorus pesticides from the environment, and counteract chemical nerve agents. It has been found that phosphotriesterase (PTE, also known as organophosphorus hydrolase, or OPH) isolated from wild-type, soil-dwelling bacteria Pseudomonas diminuta, can catalyze the hydrolysis of a broad spectrum of organophosphate triesters.6,7 PTE works particularly well for the degradation of paraoxon to diethyl phosphate and p-nitrophenol (Scheme 1) by catalyzing the hydrolysis of the P-O bond. It can also catalyze the cleavage of P-F (e.g., in sarin) and P-S bonds. Consequently, this class of natural organophosphate-degrading enzymes are of considerable interest for detoxifying organophosphate-based insecticides and chemical warfare agents.4 The remarkable enhancement of the hydrolysis of organophosphates catalyzed by the wild-type PTE can be exemplified with paraoxon as the best substrate. Raushel and coworkers determined the kinetic rate constants, kcat and kcat/Km, for the hydrolysis of paraoxon by PTE are 2280 s−1 and 6.2×107 M−1s−1, respectively, close to the diffusion-controlled limit.8 In contrast, the second order rate constant for the chemical hydrolysis of paraoxon by KOH (without PTE) is only 7.5×10−2 M−1s−1.9 Thus, at pH 7.0 the catalytic rate enhancement for the PTE-catalyzed hydrolysis of paraoxon reaches 3×1011.
One of the most significant characteristics of natural enzymes, however, is their specificity, which means each enzyme plays one or a few unique roles and works optimally only for one or a few substrates. Thus, it is not surprising at all to find that PTE is unable to tackle all unnatural phosphates and hydrolyze all the substrates at high rates. For example, while the enzyme works efficiently for the toxic paraoxon, it degrades the ubiquitous methyl parathion, chlorpyrifos, and diazinon pesticides at a very slow pace. The kinetic rate constants kcat and kcat/Km for the paraoxon analog methyl parathion, which has the phosphoryl oxygen replaced by sulfur, are only 2.4% and 0.26% of paraoxon, respectively, although the equilibrium constant Km is increased nine-fold.9 Thus, understanding the structural origin of the enzyme activity is pivotal to engineer PTE to increase its enzymatic diversity and activity.
Experimentally, both directed molecular evolution and site-directed mutagenesis techniques have been used to probe the amino acid residues responsible for the PTE catalysis and find mutated PTEs with enhanced catalytic proficiency and/or specificity.7,10–14 Cho et al.12 searched novel PTE variants by means of directed evolution and after two rounds of DNA shuffling and screening, they successfully isolated several mutants with up to 25-fold improved hydrolysis of methyl parathion. The following site-directed mutagenesis and saturation mutagenesis further confirmed the roles of the substitutions in the best mutant 22A11 where only one mutation (H257Y) is directly located in the active site.15 Griffiths and Tawfik conducted a directed evolution on PTE and found a mutant PTE-h5 which increases the rate of paraoxon hydrolysis 63 timesfaster than wild-type PTE, although the origins of the catalytic effect of PTE-h5 is still under investigation.13 By comparing the structures of PTE and acetylcholinesterase (AChE) from Torpedo california, Gopal et al. mutated residues in the binding pocket and reactive site to simulate the environment of AChE in PTE with the intention to enhance the activity of PTE towards the chemical warfare VX agent.11 Remarkably, Raushel’s group successfully increased PTE’s capability to destroy the nerve gas soman by a factor of 1000 through only three amino acid mutations (H254G, H257W, L303T).14 Recent computational simulations suggested that the enhanced efficiency of this triple mutant results from the enzyme-transition-state complementarity.16
X-ray crystallographic studies have provided critical mechanistic insights into the structure of the catalytic site and revealed that PTE is a dimeric metalloenzyme and each subunit contains a binuclear zinc center.17,18 The overall folding pattern of the monomer as shown in Fig. 1 consists of an α/β barrel with eight strands of a parallel β-pleated sheet. The two zincs, which are separated by 3.31 Å and can be substituted by other divalent metal ions such as cadmium or manganese with full catalytic activity,19–21 are situated at the C-terminal portion of the “TIM” (triosephosphate isomerase) barrel motif. Lys169, which is carboxylated, and a water molecule (or hydroxide ion) serve as bridging ligands between two zinc cations. The more buried Znα, referred to as the α-metal,21 is surrounded by His55, His57, Asp301 and Lys169, and a hydroxide ion (hydroxyl). In other words, Znα is coordinated by five ligands in a trigonal bipyramidal arrangement. The second zinc cation, Znβ, or β-metal which is more solvent-exposed,21 is tetrahedrally ligated by His201, His230, Lys169 and the bridging hydroxyl, in an approximately trigonal pyramidal arrangement. Thus, the valence of Znβ is not fully filled and there is a dangling orbital in Znβ towards the solvent-accessible surface, which remains to serve as the catalytic center in the reactions. Interestingly, the hydroxyl binds to two zinc cations in an equal distance with R(O-Zn) of 1.93 Å in the crystal structure.
Although there is an unoccupied valence bond in Znβ which will tend to form a chemical bond with a substrate, the electrostatic interaction between Znβ and the oxygen in the phosphate group in organophosphates such as paraoxon seems to be the driving force to position paraoxon into the catalytic center (Fig. 2). Surprisingly, in the crystal structure of PTE with paraoxon analog diethyl 4-methylbenzylphosphonate (EBP), Znβ and the oxygen in the phosphate group are pretty far away with the distance of 3.48 Å. But molecular dynamics simulations at the molecular mechanical (MM) level with paraoxon and sarin did show that the phosphoryl oxygen in both cases strongly coordinates with the less buried zinc ion Znβ.22 For the sarin analog diisopropyl methyl phosphonate, however, crystallographic study showed that the phosphoryl oxygen coordinates to Znβ in a much reduced distance (2.5 Å) due to the smaller molecular size compared with paraoxon.23
As for the reaction mechanism, a SN2 reaction mechanism has been postulated (Fig. 2), 18 as experiment with a chiral organophosphate substrate demonstrated that the reaction proceeds along with an inversion of stereochemical configuration at the phosphorus center.24 Crystal structures also suggested that the attacking nucleophile seems a terminally bound hydroxide, consistent with the catalytic mechanism of other binuclear metallophosphoesterases.25 Recent electron paramagnetic resonance (EPR) experiment also suggested that substrate binding via the phosphoryl oxygen at the β-metal (Mnβ in this experiment) weakens the coordination of the hydroxide bridge to the β-metal.26 Alternatively, the nucleophilic species has been proposed to be the μ-hydroxo bridge between the α- and β-metals.22,27 In general, the enzymatic reaction is postulated to be driven by the electrostatic attraction between the solvent-accessible Znβ cation and the phosphoryl oxygen. With the gradual formation of the O-Znβ bond, the hydroxyl, which initially bridges the two zinc ions, breaks its bond to Znβ and then initiates a nucleophilic attack on the phosphorus center without forming a phosphorylated enzyme intermediate. During the formation of the bond between phosphorus and hydroxide oxygen, the bond between phosphorus and p-nitrophenol oxygen breaks simultaneously. At last, paraoxon is detoxified with the release of p-nitrophenol. However, a more elaborated mechanism concerning a proton relay from the μ-hydroxo bridge nucleophile to Asp301 then to His254 and finally to Asp233 has been proposed.27 In contrast to all the above theories, Jackson et al. claimed that the μ-hydroxo bridge is not a nucleophile but acts as a base for freely exchangeable nucleophilic water molecules.28 As for the nature of the interaction between the surroundings and substrate, electrostatic theory is usually adopted to explain the significant enzymatic catalysis. But the crystal structure of PTE with substrate analogue18 reveals few electrostatic interactions between the bound inhibitor and the enzyme. Thus, the hydrophobic interactions seem the dominating force controlling the substrate selectivity and steroselectivity.29
Computationally, molecular dynamics simulations and quantum mechanical (QM) calculations have been conducted in order to elucidate the structure of PTE as well as understand the reactivity in the active site.16,22,30–36 Ko a et al.34 analyzed the large distance between the phosphoryl oxygen and the nearest zinc Znβ (3.4 Å) in the crystal structure where the substrate is the paraoxon analogue EPB, and performed molecular dynamics simulations. They found that the phosphoryl oxygen becomes strongly coordinated with Znβ if the substrate is paraoxon or sarin, in accord with the high reactivity of PTE towards these two substrates. Of particular interest to the computational chemists, however, is the nature of the bridging oxygen between two zinc atoms found in the crystal structure. Although the oxygen species was initially assumed to be a water molecule, high-level QM studies confirmed the species as a hydroxide anion.22 Similar Znα-OH-Znβ bridges exist in other metalloenzymes such as zinc β-lactamase.37 Most recently, Wong and Gao investigated the reaction mechanism of the hydrolysis of paraoxon by PTE with a dual-level combined QM/MM approach and found that the distance between two zinc ions fluctuates significantly in the reaction process.35 They further proposed a mechanism concerning no phosphodiester bridge to both zinc ions as the product intermediate which centers on the mechanism by Aubert et al..27
Due to the high computational costs, high-level QM calculations have often been performed with simplified cluster models. Kafai and Krauss determined the structure of the active site in PTE at the DFT level and the optimal cluster model is in reasonable agreement to the X-ray crystal structure.33 Krauss used a combined cluster model and effective fragment potentials (EFP) approach to determine the inherent electronic and structural characteristics of PTE,32 and explored the influence of water molecules in the active site.38 Similarly, Zheng et al. determined the structures of the active sites of Cd2+-containing PTE.30 They further performed ab initio and DFT calculations on active site models to study the reaction pathways and energy profiles of the alkaline hydrolysis of paraoxon, diisopropyl phosphorofluoridate (DFP), sarin, soman and O,O-dimethyl phosphonofluridate. Recently, Chen et al. derived the potential energy surface for the catalytic hydrolysis of paraoxon at the DFT(B3LYP) level.36
In this article, we performed extensive combined QM(PM3)/MM molecular dynamics simulations on the substrate binding, SN2 reaction, and product departure steps in the detoxification of paraoxon catalyzed by PTE, with aims to elucidate the catalytic mechanism and unveil how the well-structured micro-environment around the reactive site accelerates reactions. Single-point QM(DFT)/MM optimizations for key configurations were also conducted to examine and verify the simulation results. For comparison, simulations on the hydrolysis of diethyl p-tolyl phosphate which is derived from paraoxon by replacing the nitro group with a methyl group were also performed. The information garnered in our computations could be illuminating for the rational design of PTE mutants with improved catalytic power to not only paraoxon but other organophosphates.
PTE is a dimeric metalloenzyme and each monomer consists of 365 amino acids with about 3000 non-hydrogen atoms. Although many QM computations on small clusters of PTE with the catalytic center have been performed,22,30–34,36 the catalytic power and specificity of enzymes largely come from the overall enzyme and solvation environment.35 Thus, it is critical to consider the whole enzyme system in calculations to probe the catalytic mechanism, as high level ab initio calculations without taking the whole enzyme into account cannot offer information about the role of the bulk part of the enzyme that is not modeled.39 As a matter of fact, since most chemical reactions occur in condensed states, the impact of environmental effects and dynamics on the properties and reactivity of substances has been an intriguing research topic for computational chemists. One way to take advantage of both the accuracy of QM methods and the efficiency of MM methods is to combine QM and MM methods such that a small part of the whole system is treated explicitly by a QM method, while the remaining majority is approximated by a MM method.40 During the past decade, combined QM/MM methods represent one of the most significant advances in computational chemistry.41,42 Molecular dynamics simulations with the combined QM/MM methods can provide the details related to the catalytic activities and reaction mechanisms. For instance, although site-directed mutagenesis and kinetic isotope effect experiments afford important information on the reaction mechanism and the particular roles of residues in an active site, interpretation of these experimental data requires the assistance from theoretical studies to elucidate the energetics of enzyme catalysis. The accuracy of the combined QM/MM methods mostly relies on the following three factors: the partition of QM and MM parts for the targeted complex system, the quality of the empirical potential functions (or force fields), and the level of the QM computation. For enzyme systems, the separation of the QM and MM regions typically involves the breaking of one or more covalent bonds. In the present work, the QM part as shown in Fig. 1b includes the two divalent zinc cations and their ligands (His55, His57, Lys169, His201, His230, Asp301, and one bridging hydroxide anion), plus the substrate paraxon. An additional water molecule binding to zinc may be considered. In total, there are up to 100 QM atoms and 17019 MM atoms. Since the six amino acid residues covalently link to others, the QM/MM boundary is treated with the generalized hybrid orbital (GHO) method,43 where for each boundary atom between the QM and MM fragments a set of hybrid orbitals is dynamically determined and one of the hybrid orbitals participates in the SCF calculations for the QM region. The CHARMM-22 force field was used to represent all of the MM atoms,44 and bond lengths involving hydrogen atoms were constrained with the SHAKE algorithm.45
In terms of the QM level in combined QM/MM methods, most calculations put the QM portion at semi-empirical level (such as AM146, PM347 or SCC-DFTB42) due to the computational efficiency. Recently, density functional theory (DFT) is being adopted in combined QM/MM methods.48 This level of theory is particularly important to metalloenzymes since it is difficult to treat transition metals adequately using semi-empirical methods. In the current work, we are concerned with the hydrolysis of paraoxon catalyzed by PTE and the free energy profile for the reaction to locate the transition states as well as determine the reaction barriers. To make the calculations efficient, we used the PM3 method for the QM part. However, it is known that both the AM1 and PM3 models do not perform well for phosphorus-containing compounds, as an sp-type basis set is not enough to reproduce the properties of phosphorus compounds acceptably, and d orbitals are being added to these methods.49 Notably, the recently developed PM6 parameter set by Stewart also includes d orbitals for phosphorus.49e An alternative and cost-efficient approach is to re-optimize the parameters in comparison with high-level QM computations and derive a set of specific reaction parameters (SRP).50 Based on the reaction energy profile of paraoxon + OH− in gas phase at the B3LYP/6-31+G(d) level, we refined the PM3 phosphorus parameters (including Uss, Upp, βs, βp, Zs, Zp, and α and limited to 10% variation from the original values) and details can be found in the Supplementary Information. Since the original PM3 zinc parameters work poorly for zinc metalloenzymes,51 the updated ZnB (Zinc, Biological) parameter set developed by Brothers et al.52 was used throughout this work. We would like to reinforce that the essential value of computational simulations lies in the understanding of the enzyme catalysis, i.e., the conformational evolution as well as the accompanied energetic change along the reaction process. In our opinion, the illustration of the conformational change is more critical than the computed energy terms such as reaction barriers which heavily rely on the theoretical level used in computations and can be progressively improved by increasing the theoretical level.
As such, we chose the QM(PM3-SRP)/MM method and the software CHARMM to carry out all stochastic boundary molecular dynamics simulations.53 Crystal structure of PTE with bound substrate analog diethyl 4-methylbenzylphosphonate (PDB code: 1DPM) was used in this study as the initial configuration by replacing the inhibitor with paraoxon due to the high similarity between diethyl 4-methylbenzylphosphonate and paraoxon, though the PTE structure with diisopropyl methyl phosphonate (PDB code: 1EZ2) has a shorter distance (2.5 Å) from the phosphoryl oxygen to Znβ.23 Prior to computations, initial adjustments were made to fit paraoxon into the known pockets8 and bound state. After re-centering the middle point of two zincs and adding the hydrogen atoms using the HBUILD facility in CHARMM, we immersed the protein into a water sphere of 30 Å radius. The reaction zone is defined as a sphere of 25 Å radius, and atoms in the reaction zone are propagated according to Newtonian mechanics. Atoms within 25–30 Å comprise the buffer zone, and are retained by a harmonic restoring force and propagated by Langevin dynamics. The friction constants for protein atoms and water molecules are assigned as 200 ps−1 and 62 ps−1, respectively. Reservoir zone refers to the atoms beyond 30 Å from the center, which are fixed during the MD simulations and whose effect is simulated by deformable stochastic boundary potentials.54 For the water molecules, the TIP3P potential was used. The list of non-bonded interactions, updated every 25 steps, was truncated at 14 Å and the van der Waals and electrostatic interactions were smoothly switched off in the range of 12–13 Å. A few steps of minimizations were performed to remove the bad contacts and relax the enzyme-substrate complex. The umbrella sampling technique was adopted to generate the free energy profile or potential of mean force (PMF) by overlapping windows along a designated reaction coordinate.55 The trajectories are generated with a time step of 1 fs at constant pressure and temperature (300 K). A biasing harmonic potential with a force constant of 25~30 kcal/mol was imposed in simulations which were separated by 0.3 Å in reaction coordinates. For each simulation, the first 100–200 ps simulation brought the system to an equilibrium state, and the second 100 ps simulation generated dynamics data for further analyses.
To confirm the above QM(PM3-SRP)/MM molecular dynamics simulation results, we further performed single-point QM(DFT)/MM optimizations on selected key configurations. In these geometry optimizations, the QM region is the same as that in the above QM(PM3-SRP)/MM simulations, and the active MM atoms include all residues around two zincs within a distance of 15 Å, whereas the remaining MM atoms are frozen. The QM part was treated at the B3LYP/6-31G(d) level and the MM part was similarly described by the CHARMM force fields. An electronic embedding scheme56 incorporating the MM charges into the one-electron Hamiltonian of the QM treatment and hydrogen link atoms with charge shift model57 for the QM/MM boundary were adopted here. The ChemShell package58 integrating the Turbomole59 and DL-POLY60 programs were employed to perform these QM(DFT)/MM single-point optimizations.
The crystal structure of PTE with EBP shows that the substrate analog cannot bind to zinc effectively as the distance between the phosphoryl oxygen and the nearest zinc (Znβ) is 3.48 Å,18 so is PTE with diethyl 4-methoxyphenyl phosphate (EPO) whose distance is 3.2 Å.28 In contrast, both molecular dynamics simulations and quantum mechanical calculations confirmed that the phosphoryl oxygen can strongly coordinate with Znβ.34 Thus, we started the simulations by building a Michaelis complex where paraoxon is in the binding distance to Znβ (~2 Å). It has also been recognized that the active site of PTE is composed of three highly hydrophobic and distinct pockets, which accommodate the three substituents on the central phosphorus atom.8,18 Accordingly, we fitted the two ethoxy groups of paraoxon into the first two pockets, one of which is the large binding pocket including His257, Leu271, Phe306 and Met317, and the other is a small subunit consisting of Trp131, Ile106, Leu303, Phe306 and His57.8 The third pocket, loosely defined by Trp131, Phe132, Leu271, Phe306 and Tyr309, is suited for the leaving group, namely the nitrophenol group in our case. The side chains of these residues dictate the substrate and stereo-selectivity of PTE. With the construction of the initial Michaelis complex following the above guidelines and the addition of water sphere, we performed combined QM/MM simulations for 300 ps to generate the equilibrated conformation and used the last 100 ps to derive the average bond distances for the following discussion. Fig. 3 shows the snapshot in the end of the simulation.
While paraoxon can form a strong bond with the less buried Znβ as the O1-Znβ bond distance remains 2.03 Å with little fluctuation, the bridging hydroxyl has no sign of cleaving its bond with Znβ as the hypothesized mechanism shown in Fig. 2 suggests. As a matter of fact, the hydroxyl binds to Znα and Znβ equally in the distance of 2.00 Å. As such, both zinc ions adopt the same trigonal bipyramidal coordination pattern in the paraoxon-bound state. The Znα-Znβ distance is 3.71 Å when paraoxon is bound to Znβ, and this Znα-Znβ separation is larger than the experimental value of 3.31 Å (3.51 Å in our present study) when the substrate is not bound. In the equilibrated configuration (state B) as shown in Fig. 3, the two hydrophobic pockets 1 and 2 bind the two ethoxy groups very well via van der Waals interactions, and from our simulations we identify that the large pocket consists of Leu271, Leu272, His254, His257, Leu303 and Phe306, whereas the small pocket is composed of Gly60, Ser61, Ile106, Tyr131, Phe306 and Ser308. Gly60 was first proposed to be included in the active site by Ko a et al.,34 although its role in the selectivity of the enzyme had been confirmed experimentally.8 However, significant discrepancy between crystal structures and the present simulation results is observed for the micro-environment around the leaving group. Instead of lingering in the hydrophobic pocket 3 favored by the presumptive π-π stacking interaction, the highly polarized nitro group tends to move out and interact with solvent water molecules. As a consequence, the p-nitrophenol group partially solvates, although it does interact with a few adjacent residues such as Ala270 and Tyr309 as exhibited in Fig. 3. The polarity of the leaving group likely plays a key role in the reaction and has a direct correlation with the reaction barrier in the SN2 reaction step, as the nitrophenol part will depart in the hydrolysis step and eventually solvate in solution.
To further estimate the binding strength of paraoxon with Znβ and explain why the EBP and EPO substrates are not bound to Znβ in the crystal structure, we plotted the free energy profile for the binding process along the distance between the phosphoryl oxygen and Znβ which serves as the reaction coordinate here, i.e., RC1 = R(O1-Znβ), as shown in Fig. 4. When the distance is beyond 5.5 Å, paraoxon is essentially fully solvated and the phosphate is surrounded by water molecules to form hydrogen bonds, and the enzyme/substrate complex is in a free state (hereafter named state F). With the approach of the phosphoryl oxygen to Znβ, there are two remarkable conformational changes. One is the gradual desolvation of the phosphoryl oxygen due to the cavity size around Znβ, and the other is the evolution of the bonding pattern for Znβ which distorts from a tetrahedral coordination to a trigonal pyramidal form in preparation for the bond with the phosphoryl oxygen, as shown in Fig. 5b. Both of the conformational changes are energy-consuming which may be partially offset by the stabilizing electrostatic interaction between Znβ and O1 when they are close enough. Thus, a transition state TS1 at RC1 = 2.56 Å is observed for the binding process with a barrier of 15.2 kcal/mol. Once the substrate overcomes the barrier, there is a steep fall in energy and a sharp basin around the O1-Znβ distance of 2.03 Å is formed, indicating a strong bond between the phosphoryl oxygen and Znβ. However, the overall binding process is still endothermic by 6.8 kcal/mol, suggesting that this process is thermodynamically unfavorable. This might be the cause why most of the substrates are 3.2–3.48 Å away from the binding site in the crystal structures.
Fig. 6 illustrated the fluctuations of three key bond distances concerning the zinc ions and the bridging hydroxyl along the binding process. With the binding of paraoxon to Znβ, the distance between the two zinc ions slightly increases, from around 3.5 Å to 3.8 Å. For comparison, in the crystal structure the two metal ions are separated by 3.3 Å. When the substrate is not bound, the hydroxyl binds to both zinc ions slightly unequally and favors the less coordinated Znβ with R(Oh-Zn2) of 1.96 Å compared with R(Oh-Zn1) of 2.00 Å. In other words, the addition of one ligand to Znβ slightly weakens the Oh-Znβ bond and lengthens it from 1.96 to 2.00 Å. Interestingly, in Fig. 5a–c, we oriented the QM part except paraoxon in similar positions and it can be clearly observed that there is a rotation (counter-clockwise from the top) of the substrate to fit the two ethoxy groups into pockets 1 and 2 when paraoxon move from outside to the binding site. Similar conformational changes were detected in the simulations at the pure classical level.34
It is worthwhile to note that the recent EPR experiment on Mn/Mn-substituted PTE by Samples et al. showed that the binding of inhibitors DIMP (diisopropyl methyl phosphonate) or TEP (triethyl phosphate) to the PTE binuclear metal center diminishes the population of hydroxide-bridged species and thus weakens the coordination of the bridging hydroxide to the β-metal ion.26 But our QM(PM3-SRP)/MM simulations favor the μ-hydroxo-bridge species at the paraoxon-binding state, in accord with the recent DFT optimizations of a simplified model of the PTE active site.36 This discrepancy may result from the difference between the metal ions, as Zn2+ has a slightly smaller ionic radius (0.65 Å vs. 0.80 Å) and thus has a higher binding capability for ligands than Mn2+.
Once paraoxon is bound to Znβ, the hydroxide ion is in a distance of around 3.8 Å with the phosphorus atom (see Table 1, State B). However, the nitrophenol group is positioned in the opposite side of the hydroxyl, and thus in a perfect orientation for a typical SN2 reaction. Experimental data also indicate that hydrolysis proceeds with a net inversion of stereochemical configuration of the phosphorus center. In the SN2 reaction mechanism, the bridging hydroxyl serves as the attacking nucleophile, while the leaving group is the p-nitrophenol anion. Thus, we explored the free energy profile along the reaction path by defining the reaction coordinate (RC2) as the distance difference between the breaking bond P-O2 and forming bond Oμ-P. Fig. 7 exhibits the energy curve corresponding to the hydrolysis step of paraoxon catalyzed by PTE, which is the rate-limiting step in the overall process including binding, hydrolysis and departing. Fig. 8 plots the evolution of several key bond distances along the reaction process.
As the pre-requisite for the SN2 reaction is the approaching between the hydroxide anion and the phosphorus center whereas the former is tightly held by zinc ions, we observe a remarkable distortion of the substrate which is accompanied by a significant energy cost due to either the structural bumping between PTE and paraoxon and/or the destabilization of the substrate per se. Fig. 8 also shows that before the first transition state (TS2) at RC2 = −0.65 Å the P-O2 bond changes little and only paraoxon itself deforms and moves to approach the hydroxyl. In other words, the reduction of the Oμ-P distance dominates the reaction coordinate change, which results in other related structural variations. For instance, in the binding state the P-O1-Znβ angle is close to linear and around 170°. But in the TS2 state, it bends to about 120° and consequently the key distance between the bridging oxygen and the phosphorus atom reduces to 2.3–2.4 Å from ~3.8 Å. As the most critical outcome of this distortion, the highest reaction barrier (18.7 kcal/mol) is experienced at this stage. For comparison, DFT computations with a simplified cluster model and the continuum solvent model predicted a barrier of 10.8 kcal/mol at this step.36 However, our QM(DFT)/MM optimizations on the whole enzyme-substrate complex with the software ChemShell57 resulted in a comparable barrier of 21.0 kcal/mol (see Table 2). Fig. 9 shows the key residues surrounding paraoxon in TS2. Much similar to the binding state, at the TS2 state both ethoxy groups are still well fitted to the large and small binding pockets respectively. In the first (large) pocket, the ethoxy group has close contacts with His254, Leu271, His257 and Phe306, whereas in the second (small) pocket Tyr131, Ile106, Ser308, Gly60 and Phe306 play roles to stabilize the substitute group. For the leaving group, the nitrophenol ring rotates and has certain contact with Ser308, Tyr309, and Leu136 from the other monomer. Notably, Ser308 and Tyr309 are in the way blocking the p-nitrophenol group when paraoxon deforms to approach the bridging hydroxyl. On the other hand, the water molecules around the nitro group make the decisive contribution to stimulate the departure of the nitrophenol group. As both QM(PM3-SRP)/MM molecular dynamics simulations and QM(DFT)/MM optimizations suggest that the rate of the detoxification of paraoxon is determined by this step, the reduction of the barrier can be accomplished by either stabilizing paraoxon or reducing the repulsion between the substrate and PTE in the deformation of paraoxon. We note that Raushel’s group’s mutant targeted to soman involves His254, His257 and Leu303 in the large binding pocket.14 It would thus be expected that the mutation of Tyr309 to a much smaller residue such as glycine or serine which can interact with the nitro group may clear the way for the reorientation of paraoxon to allow the phosphorus center to approach the hydroxide ion, and thus reduce the reaction barrier and increase the enzymatic efficiency. But Raushel’s group has systematically investigated the individual role of residues in the three hydrophobic binding pockets by mutating Cys59, Gly60, Ser61, Ile106, Trp131, Phe132, His254, His257, Leu271, Leu303, Phe306, Ser308, Tyr309, and Met317 to alanine one by one.29 Unfortunately, the Y309A mutant improves little compared with the wild type (e.g., the maximum velocity for the hydrolysis of paraoxon increases modestly from 7100 s−1 to 8100 s−1, but the value of kcat/Km decreases slightly from 6.4×107 M−1s−1 to 6.1×106 M−1s−1).29 Furthermore, the mutation of Tyr309 to a phenylalanine excluded the possibility of any hydrogen bonding interaction between Tyr309 and the p-nitrophenol leaving group which was suggested by Ko a et al..34 We were also unable to spot any hydrogen bond between Tyr309 and the p-nitrophenol group in the simulations. Thus, the barrier (18.7 kcal/mol) may most likely result from the dislocation/deformation of paraoxon bound to Znβ and there is a significant loss of the binding energy between paraoxon and Znβ due to the bending of the P-O1-Znβ angle. Any hydrogen bonding from adjacent residues to O1 can potentially help stabilize the transition state TS2 and lead to the lowering of the SN2 reaction barrier. In fact, a potential hydrogen bond between NE1 (the nitrogen atom in the indole ring) of Trp131 and the phosphoryl oxygen O1 of the substrate has been identified.18,23 Directed evolution of PTE shows that this residue prevails in the more active clones.13 Our computational simulations demonstrate that at the binding state, the distance between NE1@Trp131 and O1 is around 3.3 Å, which is slightly reduced to 3.2 Å at the transition state TS2. Conformational analysis indicates that Trp131 and Phe132 are the closest residues to the phosphoryl oxygen, but NE1@Trp131 and O1 are still too far away. Thus, any mutation which can shorten the hydrogen bond distance with the phosphoryl oxygen has the potential of leading to a mutant with enhanced activity for the detoxification of organophosphates. We expect that the mutations of Trp131 to Arg, Gln, or even Tyr may be worthwhile trying. As for Phe132, similar mutations hopefully generate positive results.
Once the reaction passes the energy height, its energy curve moves downward until an intermediate state (I1) is formed at RC2 = 0.18 Å. This intermediate state is very much similar to the ion-molecule complexes for the SN2 reactions of X− + CH3X → XCH3 + X− in the gas phases.61 Significantly, at this stage the hydroxyl has completely broken its bond to Znβ (R(Znβ-Oμ) = 3.77 Å) and started to form a covalent bond with the phosphorus atom. Lack of the bridging power from the hydroxide ion, the zinc ions departs remarkably from 3.58 Å in the TS2 state to 4.49 Å in the intermediate state, which locates only 1.9 kcal/mol above the binding state. For comparison, our QM(DFT)/MM optimizations generated a value 4.3 kcal/mol. Thus, our data are notably smaller than 9.7 kcal/mol predicted based on the simplified cluster model computations.36 Another striking feature for this state is that the p-nitrophenol anion is already in the breaking point with the phosphorus center. This picture is quite similar to the PTE-substrate complex during hydrolysis captured by Jackson et al. in the crystallography analysis where the leaving group has already been in the leaving channel and the hydrolytic product diethyl phosphate (DEP) is bound in a chelating mode at Znβ (see Figs. 2c and 2d in the reference).28 As the hydrolysis goes on via an SN2-like mechanism,24 there would be a net inversion of stereochemistry at the phosphorus center with the nuclophilic attack of the hydroxide ion. Interestingly, we observed that one ethoxy group stays steadily in the small pocket, whereas the other one moves out of the large pocket to solvate itself. It is known that the substrate specificity of PTE is mainly dictated by the steric constraints within the small binding pocket.29 After a small barrier (2.2 kcal/mol) at the QM(PM3-SRP)/MM level or 8.5 kcal/mol at the QM(DFT)/MM level characterized by the transition state TS3 (see Fig. 10), the p-nitrophenol anion leaves gradually and the phosphorus regains its linearity along the O1-Znβ bond. Both zinc ions now take the tetrahedral coordination mode and their distance fluctuates around 5.3 Å (see Fig. 5f). As a matter of fact, Fig. 8 displays a good correlation between the Znα-Znβ distance and the Znβ-Oμ bond length. Interestingly, the hydroxyl group binding to phosphorus now positions itself in pocket 1, while one ethoxy group fully solvates in the solvent. The gradual decline of the energy curve highlights the high solvation energy for the p-nitrophenol anion unless it picks up a proton from a surrounding water molecule.
It should be noted that the intermediate state I1 was not identified in the recent simulations at both the QM(AM1)/MM level where both sets of parameters for P and Zn were not refined for the targeted system and the dual-level where the QM(AM1)/MM minimum energy path was corrected by the B3LYP computations of the QM part.35 But Chen et al.’s high level QM optimizations on reduced cluster models led to a similar intermediate state.36 Our result is further supported by the most recent crystallographic structure of the Co/Co-PTE-diethyl phosphate complex (PDB code: 3CAK), where the diethyl phosphate (DEP) complex symmetrically bridges the two Co2+ ions.62 This structure shows similarities with the structure I1. For comparison, in the Co/Co-PTE-DEP complex, the two metal ions are separated by a distance of 4.0 Å, and one of the phosphoryl oxygen atoms of DEP is 1.97 Å from α-metal ion, while the other oxygen is 2.17 Å from the β-metal ion. On the other hand, our I1 configuration is very similar to the transition state derived from the QM(AM1)/MM simulations, which demonstrated a very flat and broad transition plateau spanning from RC2 = −0.33 to 0.33 Å. Interestingly, the highest reaction barrier in this work (18.7 or 21.0 kcal/mol) is very close to the reaction barrier at the dual-level (18.3 kcal/mol) by Wong and Gao. For comparison, the QM(AM1)/MM level greatly overestimates the barrier (28.7 kcal/mol).35
It has been proposed that the reactivity of the μ-hydroxo bridge nucleophile could be remarkably enhanced by a proton transfer to Asp301 followed by a proton relay from Asp301 to His254 to Asp233, as the mutations to either Asp233 or His254 resulted in a decrease in the kinetic constants for paraoxon.27 However, our simulations indicate that although the hydrogen of the bridging hydroxyl is close to the carboxylate group of Asp301 and thus can form a stabilizing hydrogen bond with Asp301, the subsequent activation of the p-nitrophenol leaving group essentially requires only a little energy and no proton transfer to Asp301 is either necessary or observed in our simulations. In addition, the proton transfer from the hydroxyl to Asp301 would leave the hydrolytic product negatively charged. As a consequence, the product would bind to the binuclear metal ions closely and the subsequent departure would be difficult.
With the gradual leaving of the p-nitrophenol group, the hydrolytic product of paraoxon, namely diethyl phosphate (DEP), binds to Znβ. Apparently, DEP has a tendency to stabilize the complex by bridging both zinc ions.62 But here we are interested in the energy cost for the departure of DEP from the active site of PTE. As the simulations goes on we observed a few water molecules around the zinc ions in the range of ~4 Å and it is expected that a water molecule would eventually displace DEP. A similar phenomenon was described by Wong and Gao, who confirmed a dynamical water molecule in the average distances of 3.3 ± 0.3 and 4.5 ± 0.4 Å with Znα and Znβ, respectively.35 Thus, in the following simulations we excluded the departed p-nitrophenol group but added a water molecule to the QM part which now consists of 89 atoms in total. The reaction coordinate (RC3) in this round was defined as the distance difference between the breaking bond Znβ-O1 (the phosphoryl oxygen) and forming bond Znβ-Oμ and the initial state was labeled as P. The energy profile corresponding to the substitution of DEP by water is shown in Fig. 11.
Interestingly, with the inclusion of one water molecule into the QM region, a few steps of minimization leads to stable structures (Fig. 5g and 5i) where the water is bound to either Znα or Znβ, depending on the initial configurations. This highlights the preference of a trigonal pyramidal pentacoordination pattern over a tetrahedral coordination pattern for zinc ions.22 At this stage, there are a few viable pathways for the subsequent process, as the bound water can deprotonate first, followed by the binding of the hydroxide ion to the other zinc ion and the expulsion of DEP. Since our focus was on the hydrolysis of paraoxon, the deprotonation of water bound to zinc was not explored in this work. Instead, we investigated the energy profile regarding the replacement of DEP by a water molecule. Fig. 11 shows that the transfer of the bound water from Znα (Fig. 5g) to Znβ (Fig. 5i) must overcome a barrier of 13.5 kcal/mol. Note that the energy curve around the transition state TS4 is discrete rather than continuous. This is due to the fact that the breaking of the Znα-OH2 bond is not considered in the current reaction coordinate which is one-dimensional. A more complete and smooth energy profile can be generated by two-dimensional simulations. Thus, the present barrier 13.5 kcal/mol represents an upbound value. The water transfer is slightly endothermic by 3.4 kcal/mol. In the intermediate state I2, Znα is tetrahedrally coordinated but Znβ is in the saturated trigonal pyramidal coordination. The subsequent departure of DEP requires the breaking of the Znβ-O1 bond which is endothermic but in contrast the solvation of DEP is stabilizing. Overall the reaction barrier is 13.1 kcal/mol. We chose the state at RC3 = 2.56 Å as the pre-final state (pF, Fig. 5k) since afterwards the energy decreasing is dominated by the solvation of DEP which is highly polarizable and can strongly interact with water molecules. Considering the solvation energy of DEP, in general the replacement of DEP by water is an exothermic process. In the pF state, both zinc ions are tetrahedrally coordinated, and one water molecule binds to Znβ. Apparently, the deprotonation process can be either step-wise or concerted. In the former mechanism, the water loses one proton at Znβ, and then the hydroxide ion bends to bridge both zinc ions. In contrast, in the concerted mechanism, the deprototonation occurs simultaneously with the bridging towards Znα. In terms of the identity of the base for the deprotonation, both solvent water and Asp301 are plausible based on the stoichiometry, as the proton ultimately will neutralize the earlier p-nitrophenol anion. Wong and Gao suggested that the proton of the bound water molecule is taken up by His254 at low pH or released to the solvent at high pH.35
The major finding in the above simulations of the hydrolysis of paraoxon catalyzed by PTE is that the rate-limiting step is the preparation step for SN2 reaction where both the geometry and the binding mode of paraoxon to Znβ distort. The actual SN2 reaction proceeds fast as the barrier is only 2.2–8.5 kcal/mol. This is rationalized by the high polarity and solvation energy of the departing group, p-nitrophenol anion. To verify this mechanistic explanation, we studied the hydrolysis of diethyl p-tolyl phosphate which differs from paraoxon only by replacing the nitro group with a methyl group. Fig. 12 shows the energy profile for the replacement of the tolyl group by the bridging hydroxyl. The comparison between Fig. 12 and Fig. 7 exhibits a few interesting features. Firstly, the rate-limiting step is the same with essentially the same barriers. Structural analyses of PTE complexed with diethyl p-tolyl phosphate show that the binding state and the transition state (TS2) are much similar to Figs. 5c and 5d, indicating that barriers originate from the same source, namely the structural deformation of the active site. A prominent structural parameter is the P-O1-Znβ which evolves from nearly linear in the state B to around 120° in the TS2 state. Secondly, with reference to the binding state, the intermediate state (I1) is unstable by 8.9 kcal/mol. For paraoxon, the value is only 1.9 kcal/mol. Finally, the striking difference between paraoxon and diethyl p-tolyl phosphate lies in the nucleophilic attack by the bridging hydroxide, which is understandable due to the very different natures of the departing groups. For diethyl p-tolyl phosphate, the reaction barrier is 14.5 kcal/mol, compared with only 2.2 kcal/mol for paraoxon. This significant difference is reflected in the overall reaction rates. The nitro group is a strong electron-drawing functional group while the methyl group is a moderate electron-donating group. Since the turnover of diethyl 4-methoxyphenyl phosphate (EPO) is 2×105-fold slower than paraoxon,20,28 we expect that hydrolysis of diethyl p-tolyl phosphate would be even much slower than EPO.
Due to the variation of organophosphates, the elucidation of the catalytic mechanism of PTE on the hydrolysis of paraoxon is essential for the rational engineering of PTE to increase its efficiency and change its specialty. While a SN2-like reaction mechanism has been generally recognized,18 controversies exist in the aspect of general bases and proton transfer pathways. The present computational simulations showed that the binding of paraoxon to the solvent-accessible Znβ cation does not break the bonds of the bridging hydroxide to both zinc ions, and consequently a bridging hydroxide27 rather than a terminally bound hydroxide25 is the nucleophile to the phosphoryl center. Based on the kinetic assays and mutations of Asp233, His254 and Asp301, Aubert et al. hypothesized that the nucleophilic attack on the phosphorus center by the bridging hydroxide is simultaneously accompanied by a proton transfer from the hydroxide to Asp301, and with the departure of the leaving group, the phosphate anion bridges the two zinc ions (see Fig. 13).27 The transferred proton is shuttled away from Asp301 to His254 to Asp233 and eventually to the bulk solvent. The active site will be regenerated from a bound water molecule once the product is released or replaced by water.
The recent computational simulations by Wong and Gao,35 however, disapproved the possibility of proton transfer from the hydroxide ion to Asp301 when the former attacks the phosphoryl center. Thus, there is no intermediate state where the phosphatediester anion forms a bridge using two phosphate oxygen atoms to coordinate both zinc ions as shown in Fig. 13. Instead, they concluded that only the P=O oxygen is ligated to Znβ and the phosphodiester product is neutral in the enzyme active site.
In our present combined QM(PM3-SRP)/MM molecular dynamics computations, we re-optimized the PM3 parameters for phosphorus with reference to the DFT computations and used the updated parameters for zinc ions. The simulations of the binding step clarified that paraoxon can strongly bind to Znβ, but the overall binding process is endothermic. This explains the crystal structures where the phosphoryl oxygen binds to the zinc ion only loosely in the distance range of 3.2–3.48 Å. Similar to the Wong-Gao finding, we were unable to observe the proton transfer from the hydroxide ion to Asp301. Different from the Wong-Gao mechanism, we found that the coordination of the nucleophile hydroxide oxygen to Zn breaks rather late along the reaction path, at RC2 = −0.6 Å, as soon as the free energy shows steep increase. Of interest, we identified an intermediate where the hydroxide cleaves its bond with Znβ and binds to the phosphoryl center. This intermediate structure is very similar to the recent crystal structure of Co/Co-PTE complexed with diethyl phosphate.62 The rate-limiting step comes from the distortion of paraoxon in order to access the hydroxide oxygen, and the subsequent SN2 reaction is rapid. These findings were further confirmed by the QM(DFT)/MM optimizations on the whole enzyme-complex for key configurations. The comparative simulations of the hydrolysis of diethyl p-tolyl phosphate by PTE demonstrated that the difference lies in this SN2 reaction step. Based on the present computations, we propose a modified kinetic model27 for the catalytic mechanism by PTE on organophosphates as
where the first step refers to the reversible binding of a substrate to PTE to form the Michaelis complex (EA). Afterwards, there is a reversible preparation step for the subsequent SN2 reaction and in this step, the substrate deforms significantly in order to approach the bridging hydroxide oxygen. The nucleophilic attack to the central phosphorus atom by the bridging hydroxide ion results in the irreversible cleavage of the P-O bond and the negation of an enzyme product complex (EP). The final step is the release of the products with the regeneration of the free enzyme. Experimental data have confirmed a change in rate-limiting step depending on the pKa of the leaving group.63 For fast substrates such as paraoxon, the rate-limiting step is the preparation step, whereas for slow substrates the P-O bond cleavage is rate-liming. We note that Aubert et al. also concluded that for paraoxon the chemical step is not entirely rate-limiting, but they claimed that the rate-limiting process is the regeneration of the resting state accompanied by a proton abstraction of the nucleophilic water.27
This work was supported by the Keck foundation, the National Institute of Health (NIH) and by funds from the Faculty Research and Creative Activities Support Fund, Western Michigan University. YM also acknowledged the support from Xiamen University of China for his collaborative visits.
Supporting Information Available: Description of the PM3-SRP parameters for the phosphorus atom.