Theoretical Modeling Study for the Phosphonylation Mechanisms of

Feb 28, 2008 - Jing Wang,† Jiande Gu,*,†,‡ and Jerzy Leszczynski*,†. Computational Center for Molecular Structure and Interactions, Department...
0 downloads 0 Views 517KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 3485-3494

3485

Theoretical Modeling Study for the Phosphonylation Mechanisms of the Catalytic Triad of Acetylcholinesterase by Sarin Jing Wang,† Jiande Gu,*,†,‡ and Jerzy Leszczynski*,† Computational Center for Molecular Structure and Interactions, Department of Chemistry, Jackson State UniVersity, Jackson, Mississippi 39217, and Drug Design & DiscoVery Center, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai 201203 People’s Republic of China ReceiVed: August 30, 2007; In Final Form: NoVember 21, 2007

Potential energy surfaces for the process of phosphonylation of the catalytic triad of acetylcholinesterase by sarin have been explored at the B3LYP/6-311G(d,p) level of theory through a computational study. It is concluded that the phosphonylation process involves a critical addition-elimination mechanism. The first nucleophilic addition process is the rate-determining step. The following elimination process of the fluoride ion comprises a composite reaction that includes several steps, and it occurs rapidly by comparison with the rate-determining step. The mobility characteristics of histidine play an important role in the reaction. A double proton-transfer mechanism is proposed for the catalytic triad during the phosphonylation process of sarin on AChE. The effect of aqueous solvation has been considered via the polarizable continuum model (PCM). One concludes that the energy barriers are generally lowered in solvent, compared to the gas-phase reactions.

Introduction Acetylcholinesterase (AChE)1,2 (EC 3.1.1.7) is one of the fastest enzymes for the catalytic reaction of a serine hydrolase, which is responsible for the degradation of neurotransmitter acetylcholine (ACh). AChE has two important ligand binding sites,3-6 the peripheral site and the acylation site. The acylation site lies at the bottom of the gorge, involving a so-called catalytic triad composed of a histidine residue, a glutamic acid residue, and a serine residue (residues numbering in TcAChE read as: His440, Glu327, and Ser200). The initial step of the acylation reaction catalyzed by AChE was explored by the theoretical study at the QM/MM level. It reveals the unique role of the catalytic triad and oxyanion hole in the AChE catalysis.7 Reactions of the natural substrate and the inhibitors with AChE have been reported based on experimental and theoretical studies.8-19 Similar catalytic triad mechanisms are also demonstrated in some serine proteases.20-23 The catalytic activity of AChE may be interfered by inhibitors either through binding to the peripheral site to block the entrance of the active gorge24-27 or through targeting directly on the active site to damage the enzyme activity.28,29 Organophosphorus acid anhydrides (OPs) are essentially irreversible inhibitors of serine hydrolases.30-34 As one of the representative of wellcharacterized OP-type inhibitors, sarin (GB, i.e., O-isopropylmethylphosphonofluoridate) interacts with AChE by blocking the active center of AChE through formation of a phosphonate ester bond to the active site on serine. The phosphonylation yields a considerably more stable and complicated sarin-serine adduct, which undergoes an aging process that results in C-O bond cleavage of the phosphonate ester bond and subsequently the loss of enzyme activity of AChE.25 The X-ray structure of * To whom correspondence should be addressed. E-mail: jiandegush@ go.com (J.G.); [email protected] (J.L.). † Jackson State University. ‡ Chinese Academy of Sciences.

the aged phosphonylated AChE sheds light on the atomic-level structural explanation for the substrate selectivity of AChE.32 The stereoselectivity of the phosphonylation reaction and the effects of the adduct configuration on the aging process were examined by four stereomers of soman with two chiral centers for human AChE (HuAChE) and its selected active center mutants.35 Rate-limiting P-O fission was reported in the selfstimulated inactivation of AChE by 4-nitrophenyl 2-propyl methylphosphonate.36 The molecular dynamics (MD) simulations were used to assess the molecular origins of stereoselectivity of phosphonylation for the solution structures of the pentacoordinate and tetracoordinate PSCS and PRCS adducts of TcAChE that are formed with soman.37 By use of molecular mechanics and MD combined with semiempirical calculations, the origins and diversity of the aging reaction in phophonate adducts of serine hydrolase enzymes including enzymes of AChE, trypsin, and chymotrypsin were studied.33 Matrix-assisted laser desorption-ionization time-of-flight mass spectrometry analysis was performed to resolve aging reaction pathways for various nerve agent phosphyl adducts with the active site of AChE.38 Although phosphonylation of the active-site serine of AChE has been studied for decades, the related mechanisms still remain unclear. The biochemical mechanism of OP inactivation of AChE was suggested to be an in-line displacement process which is initiated by precursory phosphorylation at the catalytic serine residue, and the phosphorylation of AChE by an OP that occurs simultaneously with the ejection of a leaving group to yield a stable, covalent phosphoserine ester bond.33,39-42 Kovach and co-workers reported the temperature dependence of the solvent isotope effect for the catalytic recruitment in the inactivation of AChE by soman, of which the rate-limiting process is proposed with less contribution of a general-basecatalyzed heavy-atom reorganization and more contribution of an induced-fit conformational change.43 Solvolysis of sarin and a VX [O-ethyl S-(2-diisopropylamino)ethyl methyl-phospho-

10.1021/jp076974w CCC: $40.75 © 2008 American Chemical Society Published on Web 02/28/2008

3486 J. Phys. Chem. B, Vol. 112, No. 11, 2008 nothiolate] model compound were theoretically studied at the MP2/6-31+G(d)//mPW1K/MIDI! level, and the results support the initial formation of trigonal bipyramidal intermediates and demonstrate kinetic selectivity for nucleophilic attack on the opposite direction of the more apicophilic methoxide ligand.44 The reaction of chiral isomers of isoparathion methyl with solubilized rat brain cholinesterase (RBAChE) was studied to reveal the consequence of phosphorus stereochemistry upon the postinhibitory reaction kinetics of AChE poisoned by phosphorothiolates.45 Mutant cholinesterases were supposed to possess enhanced capacity for reactivation of their phosphonylated conjugates.46 The potential energy surfaces (PESs) for the phosphonylation by sarin of the active site serine of the catalytic triad of AChE have been previously studied at the B3LYP/6-311G(d,p) level of theory. The serine residue was modeled by a methanol and a methoxide ion, respectively. Such models represent the nucleophiles under two extreme reaction conditions.47 Unlike the in-line displacement mechanisms, the results (pathways A and B) suggest that the phosphonylation process favors a major addition-elimination mechanism. To further investigate how the catalytic triad activates the serine residue in the phosphonylation reaction and obtain the insights of the phosphonylation mechanisms between sarin and the catalytic triad of AChE, the more complex reaction models that include residues of histidine and glutamic acid, have been considered in the present study. The effect of aqueous solvation was accounted for via the polarizable continuum model (PCM) to evaluate the solvent influences on the phosphonylation inhibition reaction. The study might be of importance for the design of TS analogue inhibitors of pharmacologically relevant serine hydrolyses. Methods The studied models have been fully optimized by analytical gradient techniques. The density functional theory (DFT) with Becke’s three-parameter (B3)48 exchange functional along with the Lee-Yang-Parr (LYP)49,50 nonlocal correlation functional (B3LYP) was applied in this study. The DFT results are comparable with those obtained by MP2 level calculations in the studied system.47 The standard valence triple-ζ basis set, augmented by d-type polarization functions for heavy elements and p-type polarization functions for H, 6-311G(d,p),51 was used. In the analysis of harmonic vibrational frequencies, the force constants were determined analytically for all of the complexes. The stationary structures are found by ascertaining that all of the harmonic frequencies are real. The corresponding transition states were characterized by the existence of single imaginary vibrational frequencies as the saddle points on the corresponding PES. An intrinsic reaction coordinate (IRC) analysis was carried out to confirm that the transition states connect the corresponding minima. All computations were performed by using the Gaussian 03 package of programs.52 The polarizable continuum model (PCM) self-consistent reaction field of Tomasi and coworkers with a dielectric constant ∈ ) 78.39 (water)53 was applied for all gas-phase-optimized structures to evaluate the solvation effects on the reactions. In this work, the models for sarin and serine are constructed the same way as those in our prior study.47 Sarin is modeled as O-methyl methylphosphonofluridate, and the serine residue is simplified by methanol (CH3OH), which is labeled as Ser-m. The histidine residue is represented by an imidazole ring, and the glutamic acid residue is expressed by a formate anion. They are labeled as His-m and Glu-m, respectively.

Wang et al. Results and Discussions At the active site the serine residue possesses nucleophilic characteristics due to the proton transfer to the imidazole nitrogen of the histidine residue in the catalytic triad. The obtained imidazolium ion is further stabilized by the polarizing effect of the carboxylate ion of the glutamic acid residue in the triad. In our previous study, the reaction site of serine residue was modeled as Ser-m to explore the phosphonylation reaction mechanisms.47 For a better understanding of the function of the catalytic triad, it is also crucial to determine how the histidine and glutamic acid residues play their roles in the phosphonylation inhibition. Herein various reaction pathways have been explored for which modeled His-Ser and modeled Glu-His-Ser are selected as the reactants so that the function of each residue of the catalytic triad in the irreversible phosphonylation reaction can be clearly profiled. While His-Ser-m is considered in the phosphonylation reaction, two extreme reaction conditions exist, which are dipicted by pathways C and D. The C pathway describes the reaction when the proton (Hδ1) at the Nδ1 position on the imidazole ring of histidine residue has totally transferred to the glutamic acid residue. The D pathway portrays the opposite extreme situation when histidine holds the Hδ1 proton during the phosphonylation reaction. The intact catalytic triad GluHis-Ser-m is taken into account for the reaction, which is depicted as the E pathway. C Pathway. In this pathway, the histidine residue is considered to be unprotonated, which simulates the situation in which the Hδ1 proton of the imidazole ring has already been seized by the glutamic acid residue. The structures of all intermediates and transition states in the C pathway are shown in Figure 1, and the reaction coordinate is depicted in Figure 2. The energy properties of all related structures are summarized in Table 1. In the C reaction pathway, Ser-m and His-m are stabilized by hydrogen-bonding interactions and attack the phosphorus of sarin acting on the opposite side from fluoride. The Hγ proton from Ser-m is simultaneously transferred to the imidazole ring of histidine at the N2 position. The protonated His-m then easily rotates to the direction of the fluoride ion of sarin to assist with the rupture of the fluorine bond. The Hγ proton, which is transferred to histidine at the first nucleophilic step, remains bonded to the N2 of the imidazole ring until the phosphonylation reaction completes. This pathway reveals a critical additionelimination reaction. For the structure representing the energy minimum of the reactants, a strong OH‚‚‚N-type H-bond is formed between the hydrogen (Hγ) of the hydroxyl group from Ser-m and the nitrogen (N2) of the imidazole ring from histidine for the near attacking complex (NAC) INTc1, with the atomic distance of 1.60 Å and a bond angle of 174.1°. The atomic distance between P and Oγ of serine is predicted to be 3.75 Å, indicating no interactions between these two species for the NAC of INTc1. The first transition state TSc1 on the PES of C Pathway is located 9.0 kcal/mol above INTc1 (∼8 kcal/mol after zero-point correction). The single imaginary frequency of TSc1 amounts to 156i cm-1. This vibrational mode reveals the simultaneous processes of both nucleophilic attack on the oxygen (Oγ) of serine on the phosphorus atom of sarin from the backside of the fluoride ion and the proton (Hγ) transfer from Oγ of Ser-m to N2 of His-m. The atomic distance between Hγ and N2 is predicted to be 1.40 Å, implying a strong N‚‚‚HO-type H-bond. The P-Oγ atomic distance is 2.23 Å, more than 1.5 Å shorter than that of INTc1, which suggests partial bonding interaction

Phosphonylation of the Catalytic Triad of AChE

J. Phys. Chem. B, Vol. 112, No. 11, 2008 3487

Figure 1. Optimized structures of all intermediates and transition states in the C pathway at B3LYP/6-311G(d,p) level. (Atomic distances in angstroms; INT represents intermediate, and TS represents transition state; F in aqua, P in orange, O in red, C in gray and H in white balls; orange arrows illustrate the vibrational mode of the transition states.)

on this site between AChE and sarin. The bond length of P-F increases slightly (about 0.1 Å longer than that of INTc1). Meanwhile, a pseudorotation of the methyl and methoxy groups of sarin is predicted for TSc1. The Gibbs free energy of activation between INTc1 and TSc1 is calculated to be 11.4 kcal/mol (at 25 °C and 1 atm pressure). The equilibrium reaction constant of this step is therefore predicted to be 2.63 × 104 s-1. Similar to the characteristics of the A pathway,47 a relatively stable local minimum INTc2 which follows TSc1 is featured as the pentacoordinated phosphorus compound with a unique trigonal bipyramidal (TBP) structure. The P-Oγ distance is shortened to 1.75 Å, while the P-F bond length is stretched to 1.78 Å, indicating that the serine is chemically bonded to sarin as a P-Oγ bond, and the P-F bond is meanwhile weakened in the INTc2 intermediate. INTc2 is further stabilized by a strong

NH‚‚‚O-type hydrogen bond between the Hγ(N2) of the His-m and O1 of the phosphorus compound (the H-bond length is 1.61 Å). The imidazole ring of histidine is deviated slightly from the apical plane of OγPFO1. The dihedral angle of OγPO1N2 is estimated to be 30°. INTc2 is by 2.3 kcal/mol more stable than the NAC of INTc1. With almost no barrier, intermediate INTc2 transforms into intermediate INTc3 by rotation of the imidazole ring of histidine around the P-O2 bond. These intermediates are connected through the transition state TSc2 with the very low-energy barrier (0.5 kcal/mol). For INTc3, the imidazole ring of histidine is located in the apical plane of OγPFO1 (the dihedral angle of OγPO1N2 is estimated to be -179°), while the H-bond length of ((N2)Hγ‚‚‚O1) between the histidine and the phosphorus compound is stretched to 1.64 Å. INTc3 is estimated to be 5.4 kcal/mol lower in energy than INTc1. In the next step, the

3488 J. Phys. Chem. B, Vol. 112, No. 11, 2008

Figure 2. PES of the C pathway (gas phase in black and PCM model results in orange; energy units in kcal/mol).

imidazole ring of the His-m moves around the OγPFO1 plane, resulting in transition state TSc3 that is obtained after crossing the relative low-energy barrier of 2.3 kcal/mol from INTc3. The INTc4 intermediate is formed after TSc3. It is characterized by the position of C1 of histidine that is oriented to the fluoride ion. As revealed by the characteristics of vibrational mode of the single imaginary frequency of the transition state TSc4, the imidazole ring of histidine stays in the OγPFO1 plane and swings away from the phosphorus compound to reach intermediate INTc5. TSc4 lies above INTc4 and INTc5 by 2.6 and 0.7 kcal/ mol, respectively. With the assistance of both the His-m and the clip of the methyl and methoxy groups of the phosphorus compound, only 0.3 kcal/mol activation energy is needed for INTc5 to reach the transition state TSc5, whose vibrational mode of the single imaginary frequency (101i cm-1) corresponds to the rupture of the P-F bond. This is a very fast process, the corresponding equilibrium reaction constant is calculated to be 7.31 × 1012 s-1. The final reaction product INTc6 lies 5.8 kcal/mol below the NAC INTc1, which is characterized by the fluorine phosphorus distance of 3.6 Å and a strong H-bond with Hγ(N2) of histidine (with the atomic distance of 1.39 Å). Three relatively stable intermediates (INTc2, INTc3, and INTc4) in the C pathway are characterized by the TBP structures. Following the first step of nucleophilic attack, four ensuing steps were predicted. These steps however can be summarized to be a combined step for the departure of the fluoride ion, of which the highest energy barrier is predicted to be 2.6 kcal/mol. Hence the C pathway demonstrates a major addition-elimination mechanism. The first step, which is the nucleophilic attack of serine-histidine on sarin, is concluded to be the rate-determining step. Compared with the A pathway,47 which simulates the event of a complete proton transfer from the Ser-m, the C pathway has a higher activation barrier for the nucleophilic attack step than that of A pathway (9.0 vs 4.6 kcal/mol). The release of the fluoride ion becomes easier than in the corresponding step in the A pathway (2.6 vs 8.0 kcal/mol). The above discussions for the C pathway are based on the results of calculations in the gas phase. To evaluate the solvent

Wang et al. effects for this pathway, the PCM model (water is used as a solvent) was applied for all the gas-phase-optimized structures on the PES of the C pathway at the same theoretical level (B3LYP/6-311G(d,p)). As depicted in Figure 2, the PES of the C pathway with the consideration of PCM (the orange dotted line) is predicted to be similar to that of the gas phase (black line in Figure 2). The activation energy barrier for the ratedetermining step is 1.5 kcal/mol lower than that of gas phase. The rest of the reaction steps show slight differences for the energy barriers compared with those of the gas phase (less than 1 kcal/mol). It is noted that the INTc6 reaction product appears to be much more stable in polar solvent. The solvent lowers its energy below the NAC INTc1 for 12.6 kcal/mol. This suggests that the solvent to some extent may help to accelerate and facilitate the phosphonylation reaction of sarin on AChE. It should be noted that in the presence of hydrophobic groups a less polarized solvent has weaker effects on the reaction than water does. Therefore, our water and gas-phase conclusions illustrate two boundary conditions, and one expects that biological systems exhibit characteristics that are between those predicted for the two extreme cases. D Pathway. The structures of all intermediates and transition states in the D pathway are shown in Figure 3, and the PES curve is depicted in Figure 4. Corresponding to the B pathway,47 the D pathway reveals the possible phosphonylation reaction path when no proton transfer occurs from the imidazole ring of the His-m to the Glu-m of the glutamic acid residue. In the initial step of the reaction, both the His-m and Ser-m retain their own protons (Hδ1 and Hγ). The Ser-m in the D pathway first transfers its Hγ proton to histidine while it carries out nucleophilic reaction with sarin. After the His-m is deprotonated, the Hγ proton migrates and binds to O1 of sarin, and the imidazole ring of histidine turns from the nucleophilic attacking side to the fluoride ion side. The histidine is then reprotonated by Hγ. The Hγ proton later facilitates the departure of the fluoride ion from the serinesarin adduct. The NAC reactant INTd1 in the D pathway is stabilized through a hydrogen bond between N2 of the histidine imidazole ring and the Hγ hydrogen of the serine hydroxyl group. The atomic distance of N2 and Hγ(Oγ) is 1.85 Å, indicating a medium strength OH‚‚‚N hydrogen-bonding interaction. The Oγ oxygen of serine is about 3.85 Å away from the phosphorus of sarin. The activation energy barrier that needs to be overcome to reach the transition state TSd1 is predicted to be as high as 30.7 kcal/mol. TSd1 is featured by a TBP structure. Its single imaginary frequency amounts to 257i cm-1, whose vibrational mode corresponds to the nucleophilic attack of oxygen (Oγ) of serine on the phosphorus atom of sarin from the backside relative to the fluoride ion. As a consequence, the proton of the hydroxyl group of serine (Hγ) is transferred to histidine and forms a bifurcated NH‚‚‚O hydrogen bond with Oγ of serine and O1 of sarin (with the atomic distances of 1.86 Å and 1.73 Å, respectively). The imidazole ring of histidine is located in the apical plane of OγPFO1. The P-Oγ atomic distance is shortened to 1.81 Å, suggesting the bonding interaction between serine and sarin. The P-F bond is slightly weakened by elongation of its bond length to 1.73 Å. Moreover, the methyl and methoxy groups of sarin simultaneously rotate from the orientation of the attacking serine to the direction of fluoride ion. The Gibbs free energy of activation between INTd1 and TSd1 is estimated to be 35 kcal/mol (at 25 °C, 1 atm pressure). The equilibrium

Phosphonylation of the Catalytic Triad of AChE

J. Phys. Chem. B, Vol. 112, No. 11, 2008 3489

TABLE 1: Energy Properties of the Reaction Structures at the B3LYP/6-311G(d,p) Level (kcal/mol) structures INTc1a TSc1 INTc2 TSc2 INTc3 TSc3 INTc4 TSc4 INTc5 TSc5 INTc6 INTd1b TSd1 INTd2 TSd2 INTd3 TSd3 INTd4 TSd4 INTd5 TSd5 INTd6 TSd6 INTd7 INTe1c TSe1 INTe2 TSe2 INTe3 TSe3 INTe4 TSe4 INTe5 TSe5 INTe6

∆E

∆EZPE

∆G0

∆EPCM

0.0 (-1013.17900) 9.0 -11.4 0.5 -3.6 2.3 -2.2 2.5 -0.7 0.4 -5.8

C Pathway 0.0 (-1012.97998) 8.1 -0.8 -0.5 -3.8 -1.7 -3.7 -1.3 -1.9 -2.1 -8.1

0.0 (-1013.03161) 11.4 2.8 3.3 0.2 2.6 0.3 2.9 1.6 1.6 -5.4

0.0 (-1013.26182) 7.6 -1.5 -2.2 -4.0 -2.7 -4.4 -2.5 -2.9 -3.7 -12.6

0.0 (-1013.73713) 30.7 12.7 14.0 12.2 17.0 8.9 11.2 7.9 19.0 12.0 19.6 -1.4

D Pathway 0.0 (-1013.52311) 30.6 13.7 14.8 13.2 17.3 9.7 11.9 8.8 18.3 11.5 19.8 -1.5

0.0 (-1013.57411) 35 17.2 18.7 17.1 21.4 13.1 15.7 12.3 22.6 14.5 23.9 0.4

0.0 (-1013.75781) 17.2 10.9 10.5 10.3 14.7 7.8 8.6 7.6 13.7 10.2 9.1 -4.4

0.0 (-1203.04367) 13.4 9.5 10.1 8.4 8.5 5.2 5.7 1.3 4.6 -4.5

E Pathway 0.0 (-1202.80917) 11.0 10.3 10.8 8.4 6.9 6.0 6.4 2.2 5.5 -5.3

0.0 (-1202.86840) 15.2 13.8 15.1 12.9 11.7 9.7 10.4 6.4 9.8 -2.7

0.0 (-1203.12526) 12.9 9.7 10.4 8.3 8.8 6.9 7.8 4.9 7.1 -6.9

a The energy differences are calculated based on the energy value of INTc1 (in parentheses with unit of au.). b The energy differences are calculated based on the energy value of INTd1 (in parentheses with unit of au.). c The energy differences are calculated based on the energy value of INTe1 (in parentheses with unit of au.).

reaction constant of this step is calculated to be 1.22 × 10-13 s-1, indicating that this step is slow and difficult to occur. The TBP intermediate INTd2 is formed after TSd1 releases 17.9 kcal/mol of energy. The initially protonated histidine is deprotonated as the Hγ proton transfers to sarin and bonds to the oxygen O1 of sarin with the bond length of 0.99 Å. The imidazole ring of histidine shifts out of the apical plane of OγPFO1 (with the dihedral angle of N2OγPO1 to be -12°). INTd2 is further stabilized by a strong OH‚‚‚N hydrogen bond between the nitrogen (N2) of histidine and Hγ(O1) of sarin. The calculated bond distance is 1.77 Å. The bonding interactions for P-Oγ and P-F are strengthened by shorter bond lengths than those of TSd1 (1.71 vs 1.81 Å and 1.70 vs 1.73 Å, respectively). The His-m also shows mobility as the phosphonylation reaction proceeds in the D pathway. The imidazole ring of histidine turns along the P-O1 axis from the nucleophilic attacking side to the fluoride ion side. Step by step, with the hydrogen Hγ maintaining the hydrogen bonding interactions with N2 of histidine (with the atomic distances of N2‚‚‚Hγ(O1) around 1.74-1.81 Å), INTd2 prodeeds through the transition state TSd2, the intermediate INTd3, and the transition state TSd3 by relatively low-energy barriers (1.3, 1.8, 4.7, and 8.1 kcal/

mol). An intermediate INTd4 is characterized by the histidine oriented to the side of fluoride ion and the imidazole ring of histidine located in the apical plane of OγPFO1. The imidazole ring of histidine reveals the similar dynamic characteristics of the C pathway, which also moves around the plane of OγPFO1 in the D pathway. The transition state TSd4 that connects the intermediates INTd4 and INTd5 reveals that this vibrational mode has a small imaginary frequency (33i cm-1). TSd4 lies 2.3 and 3.3 kcal/mol above the INTd4 and INTd5, respectively. Intermediate INTd5 is characterized by the C1 of histidine oriented toward the fluoride ion, the same as that of INTc4 in the C pathway, and the plane of the imidazole ring of histidine again lies in the apical plane of OγPFO1Hγ. The bond length of P-Oγ is slightly shortened (1.66 Å), while the P-F bonding interaction is weakened by elongation of its bond length to 1.78 Å. This initiated the departure of the fluoride ion. INTd5 is stabilized by the H bond of N2‚‚‚Hγ(O1) with the atomic distance of 1.73 Å. From TSd1 through INTd5, all transition states and intermediates species along the D pathway are shown to have the TBP structures. By overcoming the activation energy barrier of 11.1 kcal/ mol from the INTd5, transition state TSd5 is formed with the tetrahedral structure. The His-m is re-protonated for TSd5, in

3490 J. Phys. Chem. B, Vol. 112, No. 11, 2008

Wang et al.

Figure 3. Optimized structures of all intermediates and transition states in the D pathway at the B3LYP/6-311G(d,p) level. (Symbol representations are the same as Figure 1.)

which the Hγ hydrogen moves back to N2 of histidine from O1 of sarin. The atomic distances of Hγ to O1 and N2 are predicted to be 1.41 and 1.12 Å, respectively. The bond length of P-Oγ is 1.64 Å. Meanwhile, the fluoride ion transfers away from the phosphorus (by 2.37 Å) and initiates the F‚‚‚HC hydrogen bonding interaction with the hydrogen (H1) from C1 of imidazole ring of histidine (with the atomic distance of 1.65 Å and the bond angle of FHC equal to 137°). Intermediate INTd6 is followed after the re-protonation of histidine. Two H-bonds (N2)Hγ‚‚‚O1 and (C1)H1‚‚‚F, characterized by the atomic distances of 1.73 and 1.50 Å, respectively, help to stabilize INTd6. The H1 hydrogen gets closer (by 1.05 Å) to the fluoride ion. The fluoride ion is now 3.77 Å away from the phosphorus, indicating no interactions between these two species. The imidazole ring of histidine keeps rotating counterclockwise inside the apical plane of OγPFO1, and the other transition state TSd6 is established after INTd6 overcomes

the energy barrier of 7.6 kcal/mol. The vibrational mode of the single imaginary frequency of TSd6 (206i cm-1) corresponds to the swinging of imidazole ring of histidine and the departure of fluoride ion. The departing fluoride ion forms a bifurcated hydrogen bond with Hγ and H1, for which the atomic distances are estimated to be 1.85 and 1.99 Å, respectively. The atomic distance between fluoride ion and phosphorus is 3.16 Å. Besides, another bifurcated H bond is formed for TSd6. It involves interactions between Hγ, O1, and the fluoride ion at the atomic distances of 1.92 and 1.85 Å. A stable final product INTd7 is yielded after the TSd6 releases 21 kcal/mol of energy. The H bond between the fluoride ion and H1(C1) is destroyed. The departing fluoride ion forms a HF molecule with the Hγ proton, which originates from Ser-m, implying that serine is unable to be regenerated after the phosphonylation reaction. The atomic distance between the

Phosphonylation of the Catalytic Triad of AChE

Figure 4. PES of the D pathway (gas phase in black and PCM model results in orange; energy units in kcal/mol).

fluoride ion and phosphorus is 3.78 Å. INTd7 is more stable (by1.4 kcal/mol) than the NAC INTd1. Upon inspection of the D pathway, it can be concluded that much higher energy barriers needed to be overcome than those predicted for the C pathway. The first step for nucleophilic attack is the rate-determining step. Although there are six reaction steps together, the major steps are still proposed to be two. They are the nucleophilic attack and the elimination of the fluoride ion steps, with the energy barriers of 30.7 and 11.1 kcal/mol, respectively. It is noticeable that the His-m exhibits a dynamic role in the reaction. It first turns from the orientation of nucleophilic attacking side to the fluoride ion side and then rotates its C1 to the direction of the fluoride ion. This so-called moving His mechanism54 shows significant importance in the histidine assistance for the departure of the fluoride ion. The hydrogen which facilitates the fluoride ion departure is from the initial Ser-m. The solvent effects accounted by performing the PCM model calculations on the D pathway are considerable. As shown in Figure 4, the rate-determining step with the PCM model is still the first nucleophilic attack step. Its energy barrier is sharply lowered in solvent almost to half of its value in the gas phase (17.2 vs 30.7 kcal/mol). All the other energy barriers in the D pathway are lowered by as much as 0.3-11.6 kcal/mol in solvent compared to those in gas phase. E Pathway. The entire catalytic triad Glu-His-Ser-m is taken into account as the reactant in the E pathway. Along the E pathway, six intermediates and five corresponding transition states are located, which are illustrated in Figure 5. The relative PES curve is shown in Figure 6. Along the E pathway, the Ser-m is the reacting site, and the other two residues (Glu-m and His-m) facilitate the reaction. After serine bonds to sarin on the phosphorus, Glu-m and His-m alter their orientations as His-m does in the C and D pathways. The fluoride ion subsequently departs from the serine-sarin complex with the assistance of Glu-m and His-m. A simultaneous double proton transfer is favored among the triad throughout this pathway. Both the movements of the protons among the triad and the histidine dynamic play important roles in the E pathway. The initial NAC intermediate INTe1 is formed by multiple hydrogen bonding interactions among the catalytic triad and sarin. Two strong hydrogen bonds (Glu)Oδ1‚‚‚HNδ1(His) and (His)N2‚‚‚HOγ(Ser) link the triad together, with calculated atomic distances of 1.51 and 1.70 Å, respectively. The Oγ

J. Phys. Chem. B, Vol. 112, No. 11, 2008 3491 oxygen of serine also forms a bifurcated hydrogen-bonding interaction with hydrogens from methyl and methoxy groups of sarin (with the atomic distances of 2.17 and 2.28 Å), which benefit the stabilization of the NAC. The distance between Oγ of serine and the phosphorus of sarin is predicted to be longer than 4 Å. A TBP transition state TSe1 lies 13.4 kcal/mol above the NAC INTe1. Its single imaginary frequency (646i cm-1) reveals the vibrational mode, which describes the nucleophilic attack of Oγ of serine on the phosphorus and the simultaneous double proton transfer among the catalytic triad. The Hδ1 proton is transferred from Nδ1 of His to Oδ1 of Glu, and the Hγ proton is transferred from Oγ of Ser to N2 of His. The atomic distance between Oγ(Ser) and phosphorus is 2.08 Å, and the bond length of P-F is longer (by 0.1 Å) than that of INTe1. It should be noted that the van der Waals contact distance for Oδ1(Glu) and Nδ1(His) is as short as 2.56 Å (