PROTEINS: Structure, Function, and Genetics 28:543–555 (1997) Molecular Mechanic Study of Nerve Agent O-Ethyl S-[2-(diisopropylamino)ethyl]methylphosphonothioate (VX) Bound to the Active Site of Torpedo californica Acetylcholinesterase Christine Albaret,1 Stéphane Lacoutière,1 William P. Ashman,2 Daniel Froment,1 and Pierre-Louis Fortier1* 1Centre d’Études du Bouchet, DGA, Vert le Petit, France 2U.S. Army Edgewood Research, Development and Engineering Center (ERDEC), Aberdeen Proving Ground, Maryland 21005 ABSTRACT Herein a molecular mechanic study of the interaction of a lethal chemical warfare agent, O-ethyl S-[2-(diisopropylamino)ethyl]methylphosphonothioate (also called VX), with Torpedo californica acetylcholinesterase (TcAChE) is discussed. This compound inhibits the enzyme by phosphonylating the active site serine. The chirality of the phosphorus atom induces an enantiomeric inhibitory effect resulting in an enhanced anticholinesterasic activity of the SP isomer (VXS) versus its RP counterpart (VXR). As formation of the enzymeinhibitor Michaelis complex is known to be a crucial step in the inhibitory pathway, this complex was addressed by stochastic boundary molecular dynamics and quantum mechanical calculations. For this purpose two models of interaction were analyzed: in the first, the leaving group of VX was oriented toward the anionic subsite of TcAChE, in a similar way as it has been suggested for the natural substrate acetylcholine; in the second, it was oriented toward the gorge entrance, placing the active site serine in a suitable position for a backside attack on the phosphorus atom. This last model was consistent with experimental data related to the high inhibitory effect of this compound and the difference in activity observed for the two enantiomers. Proteins 28:543–555, 1997. r 1997 Wiley-Liss, Inc. Key words: serine esterase; enantiomeric inhibition; stochastic dynamics; ab initio calculations INTRODUCTION Acetylcholinesterase (AChE) catalyzes the hydrolysis of acetylcholine (ACh) into acetate and choline at cholinergic nerve terminals.1,2 This enzyme presents a wide range of molecular forms depending on its oligomeric assembly and the way it binds to the membrane. The catalytic subunit is well conserved among the known sequences and a high r 1997 WILEY-LISS, INC. degree of identity is found for those residues implicated in the active site. Since the X-ray structure of the catalytic subunit has been solved for the Torpedo californica enzyme (TcAChE),3 it has been possible to get more insights into this region of the enzyme. The active site is found at the bottom of a deep and narrow gorge and consists of at least four domains: (1) an esteratic subsite containing the nucleophilic serine and the residues responsible for the transition state stabilization, (2) an anionic subsite that accommodates the positive choline moiety of the substrate, (3) an acyl pocket that binds the acetyl group of the substrate, and (4) a second peripheral anionic subsite (PAS), which lines the gorge entrance and is implicated in the allosteric modulation of the enzyme.4,5 Reversible bisquaternary inhibitors like decamethonium are known to derive their potency from the ability to span the two anionic subsites.6 Irreversible inhibitors generally react with the active serine to form a covalent complex.7,8 This is the case of a wide range of organophosphorus (OP) compounds that have been referenced as chemical warfare agents. They are of the type R1R2P(O)X, where X is a leaving group, and R1 and R2 are either alkyl or ester substituents. The two-step mechanism by which they inhibit AChE is described in Figure 1. Phosphonylation rate constants kp have been measured for a range of such inhibitors and have been shown to vary weakly.9–11 This is not the case of the dissociation constants Kd, which strongly depends on the nature of the groups surrounding the phosphorus atom and the absolute configuration about the phosphorus of asymmetric compounds.7,11 For this reason, most of the studies have focused on the way these agents interact with the binding pocket. We have been interested in one particularly toxic agent of this family, O-ethyl S-[2-(diisopropylamino)ethyl]methylphosphonothioate, also called VX (Fig. *Correspondence to: Dr. P.L. Fortier, Centre d’Études du Bouchet, DGA, 91710 Vert le Petit, France. Received 1 November 1996; Accepted 3 March 1997 544 C. ALBARET ET AL. Fig. 1. Two-step inhibition pathway of AChE by OP compounds. built: In the first, the leaving group of VX was positioned in the anionic subsite of TcAChE similar to that proposed for acetylcholine (model A)16,17; in the second, this leaving group was oriented toward the gorge entrance in a favorable position for a backside attack (model B). Molecular dynamics were then performed on both types of complexes, and for each enantiomer. The difference in inhibition between these enantiomers was investigated in terms of stabilization of VX in the oxyanion hole (Gly 118, Gly 119, Ala 201) and reactivity of the serine toward the phosphorus atom. For this last purpose, ab initio calculations were also performed to look over for the best orientation of the active serine Og with respect to the phosphorus atom. MATERIALS AND METHODS Fig. 2. (right). Structure of ACh (left) and VX in the SP conformation 2). This molecule displays one of the highest rate constants of inhibition,7 and the SP enantiomer (VXS) is known to be 150 times more active than its RP counterpart (VXR).12 VX has a thioic leaving group, as O-cycloheptyl S-[(trimethylamino)ethyl]methylphosphonothioate and O-isopropyl S-[(trimethylamino)ethyl]methylphosphonothioate, two related OP compounds that have been extensively studied.11,13,14 On the basis of mutagenesis experiments, it has been recently proposed that the thiocholine leaving group of these inhibitors would orientate toward the gorge entrance, allowing a backside attack on the phosphorus atom.14 This result was quite surprising, since it had been generally assumed that cationic species bound AChE in a way favoring close interactions with residues of the anionic subsite.6,15 Since no such mutagenesis studies are available for VX, we tried to ascertain if it was possible to distinguish between the two orientations by performing a molecular mechanic study of the TcAChE–VX complex. Therefore, two models were All calculations were performed on a SGI Power Indigo 2 workstation and a SGI Power Challenge 4xR8000 supercomputer. Model of TcAChE Starting coordinates of TcAChE were taken from entry 1ACE in the Brookhaven Data Bank.3,18 Hydrogens were added in INSIGHT (Biosym). All Arg, Lys, Asp, Glu, N-terminal, and C-terminal groups were assigned ionized, except Glu 443, which is located in a hydrophobic environment.3,19 Histidine protonation sites were assigned as in Gilson.19 The missing fragment 484-488 was added by using the standard Biosym library, and optimization of its geometry was achieved in DISCOVER (2.9, Biosym). Then a 200 steepest-descent steps followed by a 200 conjugategradient steps minimization was performed, constraining all heavy atoms of the crystallographic structure to their original positions (constant of 100 kcal/mol). Construction of the TcAChE–VXS Complexes VXS was built by using the program QUANTA (4.2, MSI). It was assumed to be positively charged, as 545 INTERACTION OF VX WITH ACETYLCHOLINESTERASE judged by the pKa (8.8) (D. Casagne, personal communication, 1995) of the tertiary nitrogen during enzyme inhibition. Charges were derived from ab initio calculations by using the program DMol (Biosym). Docking of the molecules was achieved in a standard way by using QUANTA (MSI) for visualization and CHARMm 23.2 for optimization. For this purpose, the phosphonyl oxygen of VX was carefully placed in the vicinity of the oxyanion hole (Gly 118, Gly 119, and Ala 201), while the phosphorus atom was positioned near the Og atom of Ser 200. The leaving group was oriented toward the anionic subsite for model A and the gorge entrance for model B. A grid on torsion angles was then performed to minimize the electronic and steric interactions with the binding pocket. Finally, a 200 steepest-descent steps, followed by a 200 adopted basis Newton Raphson steps, minimization of both complexes was achieved by keeping the protein heavy atoms fixed. The structures obtained after these minimizations were kept as initial conformations for the molecular dynamics. tions of the hydrogen bonds maintaining the catalytic triad when no hydrogen bonding term was added to the potential. Those disruptions were accompanied by a shift of the f and c angles of Ser 200 from a disallowed region of the Ramachandran plot to an authorized one. It has been shown that the unusual position of the nucleophile backbone f and c angles in the Ramachandran plot was a common feature of the a/b hydrolase fold enzymes family, and that this particular configuration around the active residue contributed to the folding of the catalytic triad.24 So we considered our modification to be pertinent for the accuracy of the molecular dynamics in our study. The van der Waals interactions between heavy atom acceptors or donors and hydrogen atoms implicated in the network were decreased such that eij for these pairs was 0.05 kcal/mol and sij was 0.4 Å,23 where eij and sij are derived from the standard Lennard-Jones combination rules: eij 5 Œeiej and sij 5 si 1 sj 2 . (1) Construction of the TcAChE–VXR Complexes VXR was built using the same procedure as VXS. For each model, the leaving group and the phosphonyl oxygen were oriented in the same direction as for the VXS. The TcAChE–VXR complexes formed resulted in a swap of the O-ethyl and methyl moieties of VXS in the respective TcAChE–VXS complexes. Minimizations of both VXR complexes were then achieved in a similar way as for the TcAChE–VXS complexes. Ab Initio Calculations After the two complexes with VXS were refined, the cartesian coordinates of the atoms of Ser 200, His 440, Glu 327, and VXS were extracted from the rest of the system. To simplify the problem, abbreviated forms of the residues were retained for the calculation20: Cb, Og, and Hg for the serine; Cb and imidazole ring for the histidine; Cg, Cd, and Oe for the glutamate. Hydrogens were added to the heavy atoms when necessary for correct hybridation. Optimization was achieved in DMol (Biosym) for the two models by using 50 minimization steps. Molecular Dynamics Stochastic boundary molecular dynamics (SBMD) were performed by using a parallelized version of the software CHARMm 23.221 purchased from MSI. The all-hydrogens model was used and water was represented by the TIP3P type.22 An explicit angledependent hydrogen bonding term was added to the potential function for the atoms of the catalytic triad implicated in the hydrogen-bond network.23 This modification was included following different molecular dynamics performed on the enzyme alone (not shown). The trajectories indeed revealed some disrup- The parameters for the complex were taken from the CHARMm parameter file except for the inhibitor for which charges were those derived from the DMol run (see above). For nonbond interactions, a cutoff of 11 Å was applied with a switching function between 8 Å and 10 Å. Dielectric constant was set to 1. The SBMD simulation method was employed as described in studies of lysozyme,25 ribonuclease A,26 and trypsin.27 Reaction region was defined by a 23 Å sphere centered on Ne2 of His440. Buffer region was made of a 2 Å layer surrounding the reaction region. The reaction zone was filled with water molecules. All nonhydrogens atoms i of the protein in the buffer region were harmonically constrained, using a boundary force defined by Fi 5 23kBTS(ri ) Dri 7Dr2i 8 (2) where kB is the Boltzmann constant, T is the temperature, Dri is the displacement vector for atom i from its position in the energy minimized structure (see above), and S(ri) is a switching function that decreases from 0.5 to 0 when going from the buffer/ reservoir boundary to the reaction region.28 The mean square displacement for atom i was calculated by using the isotropic Debye-Waller factors Bi according to the relation: Bi 5 8p2 37Dr2i 8 (3) To simplify the problem, we used the average value of B 5 19.5 Å2 calculated on the heavy atoms of the 546 C. ALBARET ET AL. Fig. 3. Stereo views of the TcAChE–VXS in model A (top) and model B (bottom). The key residues of the active site are drawn. buffer region. Solvation of the active site was achieved by using a preequilibrated sphere (25 Å radius) centered on Ne2 of His 440. Water molecules overlapping the protein–VX complex (cutoff of 2.8 Å) were deleted. Four successive rotations of this sphere were applied in order to minimize the remaining holes in the complex. For the solvent boundary forces, a radial distribution function centered on Ne2 of His 440 was applied to the oxygen atoms.28,29 The frictional coefficients related to the stochastic forces were assigned as follows: 62 ps for the heavy atoms in the protein, and 200 ps for water oxygens.29 SHAKE30 was used to constraint bonds with hydrogens, allowing a 2-fs time step. The SBMD simulation was carried out by equilibrating the system at 300 K during 10 ps. A 100-ps production dynamic followed this equilibration. Coordinates were saved every 0.4 ps. RESULTS Starting Structures The starting structures of the TcAChE–VXS complexes used for molecular dynamics are shown in Figure 3. Only the key residues of the anionic, esteric, and peripheral anionic subsites are drawn. As mentioned above, model A corresponds to a nucleophilic attack on a face adjacent to the leaving group (opposite to the methyl), and model B to a backside attack. The thioic part of VX is slightly bent in the case of model A. Distance between the phosphorus and nitrogen atoms of VX is indeed 4.29 Å for model A and 5.21 Å for model B. This difference is essentially due to a distorted conformation of the leaving group around the two carbons linking the sulfur and nitrogen atoms in model A: the S—C—C angle is 98.5°, while the C—C—N is 95.3°. Such a distorsion is not surprising, since two of the methyls of the ACh thioic group are replaced by isopropyls in the VX molecule (see Fig. 2). Both orientations show a hydrophobic interaction between one isopropyl and the ring of Phe 330. The thioic part is mainly oriented toward Trp 84 in model A, while it is close to Tyr 121 and Tyr 334, two known PAS residues,4 in model B. However, one of the two isopropyls remains in the vicinity of Trp 84 in this second model. As pointed out in the Methods section, both models were 547 INTERACTION OF VX WITH ACETYLCHOLINESTERASE TABLE I. Distances in Angstroms of O(/P)VX to N in the Oxyanion Hole for the Four Starting Complexes Complex NG118–0(/P) NG119–0(/P) N A201–0(/P) mod. A SP mod. A RP mod. B SP mod. B RP 2.41 2.52 2.57 2.47 2.74 2.75 2.71 3.00 3.23 3.21 3.75 3.90 TABLE II. Ab Initio Calculated Distances in Angstroms for the SP Conformation Model A B Fig. 4. Stereo views showing the interaction of VX with the acyl pocket in model A (top) and in model B (bottom). The inhibitor is drawn in the two conformations (VXS in black and VXR in gray). built, assuming that the phosphonyl oxygen of VX interacts tightly with the oxyanion hole through hydrogen bonds. Orientation of the leaving group was also kept for both epimers. For this reason, differences between the TcAChE–VXS and TcAChE– VXR complexes mainly arise from the positions of the methyl and O-ethyl parts of the inhibitor in the acyl binding pocket (Fig. 4). In model A, the stacking of the leaving group against the ring of Trp 84 fixes the conformation of the thioic chain in the anionic subsite for both enantiomers. Hence the O-ethyl moiety hardly accommodates the acyl pocket in the VXS conformation, as can be judged from the value of the O—C—C angle (92.3°). In the VXR conformation, however, the methyl group can bind this pocket, since its orientation in the active site is similar to that of the proposed binding site of the methyl group of the natural substrate ACh.16,17 In the case of model B, the leaving group of VX has much more room to move than in model A. As a consequence the O-ethyl part can relax more easily in the VXR conformation (the O—C—C angle is 107.2° after minimization), resulting in a greater shift of the P/O vector (Fig. 4). This shift does not modify the orientations of the isopropyl groups, which remain identical for both enantiomers. Table I summarizes Distance Ne2H440 Nd1H440 OgS200 Ne2H440 Nd1H440 OgS200 Ab initio OgS200 Oe1E327 PVX OgS200 Oe1E327 PVX 2.63 2.47 4.41 2.56 2.45 2.58 the distances of the phosphonyl oxygen to the oxyanion hole for the energetically optimized structures of the four TcAChE–VX complexes. The values are in good agreement with those recently published for the complex of TcAChE with the transition-state analogue m-(N,N,N-trimethylammonio)-2,2,2-trifluoroacetophenone (2.9, 2.9, and 3.2 Å for residues Gly 118, Gly 119, and Ala 201, respectively),17 except model B for which the hydrogen bond with Ala 201 was lost during minimization. Ab initio Calculations Geometry optimizations of the starting structures were performed in order to see what would be the optimal configuration of VX relative to the attacking serine in both models. Such gas phase calculations have to be carried out with care, since phenomena as important as solvation of the active site is not taken into account in the procedure.31 For our purpose, we reduced the system to VX and the catalytic triad (Ser 200, His 440, and Glu 327), as we were only interested by the orientation of OgS200 with respect to the phosphorus atom. Optimizations were run for VXS in the two orientations, in order to reduce the calculation time. Similar results might be expected for VXR, since the only major change is the swap of the O-ethyl and methyl moieties. In the case of model A, the distance from OgS200 to the phosphorus was 4.41 Å, while it was 2.58 Å in model B (Table II). The greater value for model A is due to the specific configuration of VX in this orientation: the sulfur atom is indeed close to OgS200 (2.65 Å) in the starting system, and optimization resulted in a shift of both sulfur and phosphorus atoms. Apart from this shift, no significant differences in the global orientation of OgS200 with respect to the phosphorus atom, as described in the last section, were observed. Concerning the hydrogen bonds implicated in the catalytic triad, values near 2.60 Å and 2.46 Å were, respec- 548 C. ALBARET ET AL. Fig. 5. Isotropic rms fluctuations for residues in the active site of TcAChE (AS, anionic subsite; PAS, peripheral anionic subsite; AP, acyl pocket; ES, esteric subsite). The results are given for the four SBMD simulations. tively, obtained for the OgS200–NeH440 and the Nd1H440–Oe1E327 distances in the two optimizations (Table II). The latest value is close to the one of the crystal structure (2.5 Å). Most surprising is the value for the OgS200–NeH440 distance that is significantly lower than the 3.1 Å value reported by Susmann and colleagues.3 It is most likely that binding of the inhibitor induces a compression of the catalytic triad allowing the proton transfer between Ser 200 and His 440 to occur.21,32 Molecular Dynamics A molecular dynamic study was performed in order to probe the behavior of the four complexes. The quality of the SBMD simulation was assessed by comparison of the root-mean-square (rms) fluctuations of the main residues implicated in the active site from the dynamics with the fluctuations from the crystal structure. The results for the four simulations are shown in Figure 5. Though the rms fluctuations are a bit larger in the dynamics, they remain close to those of the crystal structure. The largest differences affect the trajectory of the TcAChE–VXS complex in model A for which residues of the este- ratic subsite (ES) present greater fluctuations. This is due to the fact that the O-ethyl part of VX hardly accommodated the acyl pocket during the simulation and thus moved in the vicinity of the active serine, inducing larger motions of residues Gly 118, Gly 119, Ser 200, Ala 201, and His 440. It must be noted that in the four simulations residues belonging to this pocket presented few fluctuations. This suggests that the acyl pocket is minimally flexible, the Phe 288, Phe 290, Trp 233, and Phe 331 aromatic rings restricting the size of the substituents able to accommodate it. By contrast, residues of the anionic subsite were found more mobile. In particular, Trp 84 and Phe 330 presented greater fluctuations in simulations of model B. Trp 84 belongs to the V loop of TcAChE,3 which is implicated in the allosteric modulation of the enzyme.33 Thus, the mobility of this residue is not surprising. Phe 330 ring rotated to stack against the leaving group of VX. Residues of the peripheral anionic subsite (PAS) were less perturbated in the simulations of model B than in the simulations of model A. For instance, hydrogen bond between the hydroxyl hydrogen of Tyr 334 and Od1 of Asp 72 was well conserved during both dynamics of 549 INTERACTION OF VX WITH ACETYLCHOLINESTERASE Fig. 6. Time series from the four trajectories for the distances between OgS200 and Ne2H440 (top) and between Nd1H440 and Oe1E327 (bottom). Model A is depicted in gray, model B in black. TABLE III. Average Values and rms Fluctuations in Angstroms for the Distances in the Catalytic Triad 7r8 Model A B Distance Ne2H440 Nd1H440 Ne2H440 Nd1H440 OgS200 Oe1E327 OgS200 Oe1E327 rms rmax rmin SP RP SP RP SP RP SP RP 2.94 2.80 2.55 2.62 3.03 2.93 3.37 2.62 0.11 0.12 0.07 0.16 0.15 0.09 0.07 0.14 2.69 2.56 2.39 2.46 2.66 2.72 2.99 2.38 3.29 3.13 2.87 2.81 3.41 3.17 3.84 3.02 model B, although these two residues remained close to the leaving group of VX. Deviation of Tyr 121 in the TcAChE–VXR complex dynamic in model A is due to the proximity of the O-ethyl moiety of VX in this particular conformation. In all trajectories, the catalytic triad was nearly conserved (see Fig. 6 and Table III). This suggests that the modifications we introduced in the potential function of CHARMm in order to parametrize the catalytic triad hydrogen bonds were relevant. The greater fluctuations were observed for the OgS200–NeH440 distance in the trajectory of the TcAChE–VXR complex in model B. As for the simulation of the TcAChE–VXS complex in model A, this perturbation is due to the motion of the O-ethyl part of VX. The specific mode of binding of VX in model B makes that this moiety in VXR makes close van der Waals contacts with the OgS200 atom. This resulted in a weakening of the OgS200—NeH440 hydrogen bond during the simulation. We looked forward to the interaction of the phosphonyl oxygen of VX with the oxyanion hole (Gly 118, Gly 119, Ala 201). The corresponding hydrogen bonds network has been indeed reported to account for more than 30% of the free energy in the transition-state stabilization of the TcAChE–TMFPA complex.34 In all simulations the interaction with Ala 201 was almost absent (see Fig. 7 and Table IV), suggesting a minor role for this residue in the 550 C. ALBARET ET AL. Fig. 7. Time series from the four trajectories for the distances between N of Gly 118 (top), Gly 119 (middle), Ala 201 (bottom), and the phosphonyl oxygen of VX. Model A is depicted in gray, model B in black. recognition process of VX. The TcAChE–VXS complex dynamic in model B was particularly satisfactory since the two other hydrogen bonds with NG118 and NG119 were conserved. This was not the case in the three remaining trajectories for which at least one of these two hydrogen bonds was lost during the simulation. For the TcAChE–VXR complex in model B, NG118 remained far from the phosphonyl oxygen while NG119 interacted closely with it. In the case of model A, hydrogen bonds with NG118 and with NG119 were, respectively, broken in the dynamics of the TcAChE–VXS and the TcAChE–VXR complexes. Each of these hydrogen bonds remained weak. In conclusion, VXR in model B was less stabilized in the oxyanion hole than VXS, whereas similar behaviours were observed for both conformations in model A. Concerning the nucleophilic attack, no significant modifications in the global orientation of each enantiomers were observed during the four simulations (Fig. 8 and Table V). The value of the OgS200-P–S 551 INTERACTION OF VX WITH ACETYLCHOLINESTERASE TABLE IV. Average Values and rms Fluctuations in Angstroms for the Distances in the Oxyanion Hole 7r8 Model A O(/P)VX B O(/P)VX rms rmin rmax Distance SP RP SP RP SP RP SP RP NG118 NG119 NA201 NG118 NG119 NA201 3.38 2.75 3.84 3.09 2.56 4.02 3.08 3.60 3.66 4.04 2.74 3.84 0.20 0.10 0.10 0.10 0.06 0.13 0.09 0.13 0.10 0.17 0.10 0.13 2.73 2.53 3.55 2.97 2.42 3.71 2.81 3.22 3.36 3.64 2.51 3.59 3.82 3.07 4.11 3.50 2.75 4.45 3.32 3.91 3.90 4.50 3.03 4.11 Fig. 8. Time series from the four trajectories for the distance between OgS200 and the phosphorus atom of VX (top) and for the angle of the nucleophilic attack (bottom). Model A is depicted in gray, model B in black. angle in the two dynamics of model A, around 50°, may appear curious. In fact, the leaving group of VX is perpendicular to the attacking nucleophile, but the distorded conformation of the chain makes that the OgS200-P–S angle is 30° lower than the 90° expected value for an adjacent attack. In the case of model B, Ser 200 remained in a favorable position for a backside attack for both enantiomers (OgS200-P– S < 170°). However, the attacking oxygen was ,3Å farther from the phosphorus atom of VX in the TcAChE–VXR complex simulation, the value for this distance in the TcAChE–VXS complex being very close to the ab initio calculation result. DISCUSSION The chiral nature of the complexes of asymmetric OP compounds with AChE has long been a subject of interest (see the extensive review of de Jong et al.7). The great differences observed in the inhibition constants and toxicity for the two enantiomers of a 552 C. ALBARET ET AL. TABLE V. Average Values and rms Fluctuations in Angstroms for the Distance Between the Attacking Oxygen and the Phosphorus Atom of VX 7r8 Model A B Distance OgS200 OgS200 PVX PVX rms rmax rmin SP RP SP RP SP RP SP RP 4.86 3.02 3.90 4.84 0.11 0.11 0.14 0.11 4.62 2.74 3.42 4.54 5.18 3.34 4.39 5.09 same OP compound seems to be related to structural determinants of the AChE active site at the binding stage.11,14,35 There are indeed little changes in the phosphonylation rate constants for the compounds for which they have been measured.9–11 The elucidation of the X-ray structure of Torpedo californica AChE (TcAChE)3 has allowed investigations of structure–activity relationship by molecular modeling.4,15,32,33,35–38 This technique can give useful insights in the inhibition pathway of this type of compounds, since the Michaelis complex cannot be observed due to the fastness of the reaction of phosphonylation. We have been interested in the inhibitory properties of O-ethyl S-[2(diisopropylamino)ethyl]methylphosphonothioate, also called VX. This asymmetric OP compound belongs to a family containing the most powerful chemical warfare agents ever synthesized.12 The tetrahedral arrangement about the phosphorus atom makes it a good candidate for a transition state analogue. Since the structure of a complex of TcAChE with such a transition state analogue (m-(N,N,N-trimethylammonio)-2,2,2-trifluoroacetophenone, also called TMFPA), has recently been solved,17 some assumptions concerning the structural elements of importance in the catalysis have been confirmed. Among them the existence of a tripartite hydrogen-bond network including Gly 118, Gly 119, and Ala 201, responsible for the stabilization of the carbonyl oxygen in the transition state has been assessed. It has also been emphasized that the quaternary ammonium of this compound could be stabilized by electrostatic interactions with p electrons of some aromatic residues belonging to the anionic subsite among which Trp 84 seems to play a crucial role. Such interactions had been predicted before the elucidation of the crystallographic structure.39,40 It is also known that replacement of the equivalent tryptophan in the human enzyme (HuAChE), Trp 86, by an alanine significantly alters its reactivity toward a variety of charged substrates or inhibitors while there is little effect on their noncharged homologues.15 Our study of the complex of TcAChE with VX suggests that this agent binds to the enzyme in a different way from these compounds. In placing the bulky charged leaving group in the anionic subsite (model A), no major difference could be found in the simulations of the TcAChE–VXS and TcAChE–VXR complexes allowing to explain the difference of activ- ity of the two epimers. On the other hand, when placing the leaving group of VX in the gorge entrance, the enantiomers exhibited a quite different behavior with respect to the oxyanion hole. In the case of VXS, the phosphonyl oxygen was found to make close hydrogen bonds with Gly 118 and Gly 119 during the 100 ps simulation, Ala 201 being less implicated in the stabilization of the complex. By contrast, simulation of the TcAChE–VXR complex revealed that the O-ethyl moiety of VX could not accommodate the acyl pocket and thus pushed the P/O apart from the oxyanion hole. This resulted in the break of the hydrogen bond with Gly 118, inducing most probably the lesser activity of VXR in inhibiting AChE. For this reason, we believe this model of binding (model B) more consistent than the first one (model A). Another evidence argues for model B: the leaving group is indeed oriented in a suitable direction for a backside attack. This is the most favorable case for such a nucleophilic displacement since the reaction can proceed by a direct in-line mechanism. Ab initio calculations showed that binding of VXS in this orientation resulted in a compression of the catalytic triad,20,32 while the nucleophilic oxygen and the phosphorus were found to come closer. This compression was also observed in model A, but to a lesser extent (see Table II). However, the nucleophilic oxygen and the phosphorus remained far from each other in this model, making the reaction difficult to achieve. Moreover, the nucleophile should enter the tetrahedron at a face adjacent to the leaving group, a very unfavorable case, since the pentacoordinate intermediate would have to undergo at least one Berry pseudorotation41 (or turnstile rotation42) before the release of the leaving group.43 Though the intrinsic energy barrier for this isomerisation is relatively low (5–10 kcal/mol), this represents a supplementary energy cost to be provided by the enzyme during the reaction.44 We can also note that in the case of enzymes catalyzing phosphoryl transfers, such a mechanism has never been evidenced.45 The unusual orientation of the cationic leaving group of VX in the active site should not surprise. It has been indeed already evidenced for two other compounds of the alkyl methylphosphonothioates family (O-cycloheptyl and O-isopropyl S-[(trimethylamino)ethyl]methylphosphonothioate) using mutagenesis of mouse AChE residues belonging to the INTERACTION OF VX WITH ACETYLCHOLINESTERASE 553 Fig. 10. Model of interaction of methylphosphonothioates with AChE. The case of the SP configuration is depicted. The O-alkyl and methyl moieties must be swapped to obtain the RP configuration. Fig. 9. Closeup of the active site of the TcAChE–VXS complex average structure in model B. Residues within a distance of 3.5 Å from any atom of VX are drawn. acyl pocket.14 More recently, Asp 74 (72 in TcAChE), a key residue of the peripheral anionic subsite of mouse AChE, has been shown to be responsible for the high reactivity of these two compounds. This residue interacts with the quaternary ammonium of both inhibitors through coulombic forces.44 These forces represent the main contribution in the electrostatic potential, since the orientation of the quaternary ammonium is not adequate for optimizing the cation-p interactions with the indole ring of Trp 86 (84 in TcAChE).44,46 Nevertheless, as the O-alkyl moieties of both compounds point toward the anionic subsite in the SP configuration, this tryptophan can interact with these substituents through London dispersion forces. This alternate mode of binding has been highlightened by Nair and colleagues47 using quantitative structure–activity relationship (QSAR) calculations on a series of trifluoroketone transitionstate analogues bound to the active site of TcAChE. In our study, Asp 72 was found to remain relatively close to the nitrogen atom of VX in both simulations of model B. Average values of 6.07 Å and 6.68 Å were, respectively, obtained for the Oe1D72–NVX distance in the TcAChE–VXS and TcAChE–VXR dynamics, allowing a substantial electrostatic contribution of Asp 72 to the binding of the inhibitor. Role of the anionic subsite aromatic residues in the recognition process might not be depreciated otherwise. Indeed, Trp 84 appears to remain close to the isopropyls of VX during both dynamics, while Phe 330 was found to stack against the leaving group. These two close interactions are depicted for the TcAChE–VXS complex average structure on the MOLSCRIPT plot48 of Figure 9. Thus it is most likely that this particular compound enhances its activity by optimizing London and cation-p interactions of its leaving group with the aromatic residues of the anionic locus, and by interacting electrostatically with Asp 72. CONCLUSION Using molecular dynamics together with ab initio calculations, we were able to propose a model of interaction of O-ethyl S-[2-(diisopropylamino)ethyl]methylphosphonothioate (VX) with TcAChE consistent with the high inhibitory effect of this compound and the difference of activity observed for its epimers. The leaving group points toward the gorge entrance, favoring a backside attack of the active serine on the phosphorus atom. This orientation has already been evidenced for two other O-alkylmethylphosphonothioates compounds.14 Therefore, it is most likely that inhibitors related to VX adopt the orientation depicted in Figure 10 when binding to the active site of the enzyme. Since steric limitations of the acyl pocket may prevent the accommodating of too bulky O-alkyl moieties in the RP conformation, this model is probably valid for a narrow range of compounds for which the enantiomeric effect is significant. ACKNOWLEDGMENT We thank Dr. F. Bontems for helpful discussions related to molecular dynamics. REFERENCES 1. Rosenberry, T.L. Acetylcholineserase. Adv. Enzymol. 43:103– 218, 1975. 554 C. ALBARET ET AL. 2. Taylor, P. ‘‘Pharmacological Basis of Therapeutics.’’ New York: Macmillan, 1985:110–129. 3. Sussman, J.L., Harel, M., Frolow, F., Oefner, C., Goldman, A., Toker, L., Silman, I. Atomic structure of acetylcholinesterase from Torpedo californica: A prototypic acetylcholinebinding protein. Science 253:872–879, 1991. 4. Barak, D., Kronman, C., Ordentlich, A., Ariel, N., Bromberg, A., Marcus, D., Lazar, A., Velan, B., Shafferman, A. Acetylcholinesterase peripheral anionic site degeneracy conferred by amino acid arrays sharing a common core. J. Biol. Chem. 264:6296–6305, 1994. 5. Barak, D., Ordentlich, A., Bromberg, A., Kronman, C., Marcus, D., Lazar, A., Ariel, N., Velan, B., Shafferman, A. Allosteric modulation of acetylcholinesterase activity by peripheral ligands involves a conformational transition of the anionic subsite. Biochemistry 34:15444–15452, 1995. 6. Harel, M., Schalk, I., Ehret-Sabatier, L., Bouet, F., Goeldner, M., Hirth, C., Axelsen, P.H., Silman, I., Sussman, J.L. Quaternary ligand binding to aromatic residues in the active-site gorge of acetylcholinesterase. Proc. Natl. Acad. Sci. U.S.A. 90:9031–9035, 1993. 7. de Jong, L.P.A., Benschop, H.P. Biochemical and toxicological implications of chirality in anticholinesterase organophosphates. In ‘‘Stereoselectivity of Pesticides: Biological and Chemical Problems.’’ Ariëns, E.J., van Rensen, J.J.S., Welling, W. (eds.). Amsterdam: Elsevier, 1988:109–147. 8. Nair, H.K., Lee, K., Quinn, D.M. m-(N,N,N-Trimethylammonio)trifluoroacetophenone: A femtomolar inhibitor of acetylcholinesterase. J. Am. Chem. Soc. 115:9939–9941, 1993. 9. de Jong, L.P.A., van Dijk, C. Inhibition of acetylcholinesterase by the enantiomers of isopropyl S-2-trimethylammonioethyl methylphosphonothiate iodide: Affinity and phosphorylation constants. Biochem. Biophys. Acta 268:680–689, 1972. 10. Wustner, D.A., Fukoto, T.R. Affinity and phosphonylation constants for the inhibition of cholinesterases by the optical isomers of O-2-butyl S-2-(dimethylammonium)ethyl ethylphosphonothioate hydrogen oxalate. Pestic. Biochem. Physiol. 4:365–376, 1974. 11. Berman, H.A., Leonard, K. Chiral reactions of acetylcholinesterase probed with enantiomeric methylphosphonothioates. J. Biol. Chem. 264:3942–3950, 1989. 12. Benschop, H.P., de Jong, L.P.A. Nerve agent stereoisomers: Analysis, isolation, and toxicology. Acct. Chem. Res. 21:368– 374, 1988. 13. Berman, H.A., Decker, M.M. Chiral nature of covalent methylphosphonyl conjugates of acetylcholinesterase. J. Biol. Chem. 264:3951–3956, 1989. 14. Hosea, N.A., Berman, H.A., Taylor, P. Specificity and orientation of trigonal carboxyl esters and tetrahedral alkylphosphonyl esters in cholinesterase. Biochemistry 34:11528–11536, 1995. 15. Ordentlich, A., Barak, D., Kronman, C., Flashner, Y., Leitner, M., Segall, Y., Ariel, N., Cohen, S., Velan, B., Shafferman, A. Dissection of the human acetylcholinesterase active center determinants of substrate specificity. J. Biol. Chem. 268:17083–17095, 1993. 16. Harel, M., Kleywegt, G.J., Ravelli, R.B.G., Silman, I., Sussman, J.L. Crystal structure of an acetylcholinesterasefasciculin complex: Interaction of a three-fingered toxin from snake venom with its target. Structure 3:1355–1366, 1995. 17. Harel, M., Quinn, D.M., Nair, H.K., Silman, I., Sussman, J.L. The X-ray structure of a transition state analog complex reveals the molecular origins of the catalytic power and substrate specificity of acetylcholinesterase. J. Am. Chem. Soc. 118:2340–2346, 1996. 18. Berstein, F.C., Koetzle, T.F., Williams, G.J.B., Meyer, E.F. Jr., Brice, M.D., Rogers, J.R., Kennard, O., Shimanouchi, T., Tasumi, M. The Protein Data Bank: A computer-based archival file for macromolecular structures. J. Mol. Biol. 112:535–542, 1977. 19. Gilson, M.K., Straastsma, T.P., McCammon, J.A., Ripoll, D.R., Faerman, C.H., Axelsen, P.H., Silman, I., Sussman, 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. J.L. Open ‘‘back door’’ in a molecular dynamics simulation of acetylcholinesterase. Science 263:1276–1278, 1994. Daggett, V., Schröder, S., Koolman, P. Catalytic pathway of serine proteases: Classical and quantum mechanical calculations. J. Am. Chem. Soc. 113:8926–8935, 1991. Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swaminathan, S., Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comp. Chem. 4:187–217, 1983. Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., Impey, R.W., Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935, 1983. Brooks III, C.L., Karplus, M. Solvent effects on protein motion and protein effects on solvent motion. J. Mol. Biol. 208:159–181, 1989. Ollis, D.L., Cheah, E., Cygler, M., Dijkstra, B., Frolow, F., Franken, S.M., Harel, M., Remington, S.J., Silman, I., Schrag, J., Sussman, J.L., Verschueren, K.H.G., Goldman, A. The a/b hydrolase fold. Protein Eng. 5:197–211, 1992. Brooks III, C.L., Brünger, A., Karplus, M. Active site dynamics in protein molecules: A stochastic boundary molecular-dynamics approach. Biopolymers 24:843–865, 1985. Brünger, A.T., Brooks III, C.L., Karplus, M. Active site dynamics of ribonuclease. Proc. Natl. Acad. Sci. U.S.A. 82:8458–8462, 1985. Nakagawa, S., Yu, H.-A., Karplus, M., Umeyama, H. Active site dynamics of acyl-chymotrypsin. Proteins 16:172–194, 1993. Brooks III, C.L., Karplus, M., Pettitt, B.M. ‘‘Proteins. A Theoritical Perspective of Dynamics, Structure, and Thermodynamics.’’ New York: Wiley-Interscience, 1988. Brünger, A.T., Brooks III, C.L., Karplus, M. Stochastic boundary conditions for molecular dynamic simulations of ST2 water. Chem. Phys. Lett. 105:495–500, 1984. Ryckaert, J.-P., Ciccotti, G., Berendsen, H.J.C. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. comput. Phys. 23:327–341, 1977. Mulholland, A.J., Grant, G.H., Richards, W.G. Computer modelling of enzyme catalysed reaction mechanisms. Protein Eng. 6:133–147, 1993. Benscura, A., Enyedy, I.Y., Kovach, I.M. Probing the active site of acetylcholinesterase by molecular dynamics of its phosphonate ester adducts. J. Am. Chem. Soc. 118:8531– 8541, 1996. Ordentlich, A., Barak, D., Kronman, C., Ariel, N., Segall, Y., Velan, B., Shaffermann, A. Contribution of aromatic moities of tyrosine 133 and of the anionic subsite tryptophan 86 to catalytic efficiency and allosteric modulation of acetylcholinesterase. J. Biol. Chem. 270:2082–2091, 1995. Sussman, J.L., Harel, M., Silman, I. Three-dimensional structure of acetylcholinesterase and of its complexes with anticholinesterase drugs. Chem. Biol. Interact. 87:187– 197, 1993. Ordentlich, A., Barak, D., Kronman, C., Ariel, N., Segall, Y., Velan, B., Shafferman, A. The architecture of human acetylcholinesterase active center probed by interactions with selected organophosphate inhibitors. J. Biol. Chem. 271:11953–11962, 1996. Quian, N., Kovach, I.M. Key active site residues in the inhibition of acetylcholinesterase by soman. FEBS Lett. 336:263–266, 1993. Radic, Z., Pickering, N.A., Vellom, D.C., Camp, S., Taylor, P. Three distinct domains in the cholinesterase module confer selectivity for acetyl- and butyrylcholinesterase inhibitors. Biochemistry 32:12074–12084, 1993. Bencsura, A., Enyedy, I., Kovach, I.M. Origins and diversity of the aging reaction in phosphonate adducts of serine hydrolase enzymes: What characteristics of the active site do they probe? Biochemistry 34:8989–8999, 1995. Höltje, H.-D., Kier, L.B. Nature of anionic or a-site of cholinesterase. J. Pharm. Sci. 64:418–420, 1975. Dougherty, D.A., Stauffer, D.A. Acetylcholine binding by a INTERACTION OF VX WITH ACETYLCHOLINESTERASE 41. 42. 43. 44. synthetic receptor: implications for biological recognition. Science 250:1558–1560, 1990. Berry, R.S. Correlations of rates of intramolecular tunneling processes, with application of some group V compounds. J. Chem. Phys. 32:933–938, 1960. Ugi, I., Marquarding, D., Klusacek, H., Gillespie, P., Ramirez, F. Berry pseudorotation and turnstile rotation. Acct. Chem. Res. 4:288–296, 1971. Gillespie, P., Ramirez, F., Ugi, I., Marquarding, D. Displacement reactions of phosphorus (V) compounds and their pentacoordinate intermediates. Angew. Chem. Int. Ed. Engl. 12:91–119, 1973. Hosea, N.A., Radic, Z., Tsigelny, I., Berman, H.A., Quinn, D.M., Taylor, P. Aspartate 74 as a primary determinant in 45. 46. 47. 48. 555 acetylcholinesterase governing specificity to cationic organophosphonates. Biochemistry 35:10995–11004, 1996. Fersht, A. ‘‘Enzyme Structure and Mechanism.’’ New York: W.H. Freeman, 1985:221–247. Dougherty, D.A. Cation-p interactions in chemistry and biology: A new view of benzene, Phe, Tyr, and Trp. Science 271:163–168, 1996. Nair, H.K., Seravalli, J., Arbuckle, T., Quinn, D.M. Molecular recognition in acetylcholinesterase catalysis: Freeenergy correlations for substrate turnover and inhibition by trifluoro ketone transition-state analogs. Biochemistry 33:8566–8576, 1994. Kraulis, P.J. MOLSCRIPT: a program to produce both detailed and schematic plots of protein structure. J. Appl. Crystallogr. 24:946–950, 1991.