How is acetylcholinesterase phosphonylated by soman? An ab initio QM/MM molecular dynamics study.

Acetylcholinesterase (AChE) is a crucial enzyme in the cholinergic nerve system that hydrolyzes acetylcholine (ACh) and terminates synaptic signals by reducing the effective concentration of ACh in the synaptic clefts. Organophosphate compounds irreversibly inhibit AChEs, leading to irreparable damage to nerve cells. By employing Born-Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, a state-of-the-art approach to simulate enzyme reactions, we have characterized the covalent inhibition mechanism between AChE and the nerve toxin soman and determined its free energy profile for the first time. Our results indicate that phosphonylation of the catalytic serine by soman employs an addition-elimination mechanism, which is highly associative and stepwise: in the initial addition step, which is also rate-limiting, His440 acts as a general base to facilitate the nucleophilic attack of Ser200 on the soman's phosphorus atom to form a trigonal bipyrimidal pentacovalent intermediate; in the subsequent elimination step, Try121 of the catalytic gorge stabilizes the leaving fluorine atom prior to its dissociation from the active site. Together with our previous characterization of the aging mechanism of soman inhibited AChE, our simulations have revealed detailed molecular mechanistic insights into the damaging function of the nerve agent soman.

Introduction

The physiological role of acetylcholinestrase (AChE) is to effectively terminate cholinergic signals by hydrolyzing acetylcholine (ACh).[1,2] AChE has a general α/β hydrolase globular fold, as illustrated in Figure 1a. Its catalytic site, which consists of a catalytic triad (Ser200, His440, and Glu327) and a three-pronged oxyanion hole (peptidic NH groups of Gly118, Gly119, and Ala201), sits at the bottom of a 20 Å deep, narrow, and highly electronegative catalytic gorge (Figure 1b).[3] (Torpedo californica AChE numbering is used throughout this paper.) Residues lining the catalytic gorge help position and stabilize the substrate within the active site for hydrolysis. Specifically, the catalytic anionic subsite (Trp84, Phe330, and Glu199) stabilizes the positively charged quaternary trimethylammonium tailgroup of ACh by cation−π and electrostatic interactions,[3,4,23,26] while the acyl pocket (Trp233, Phe228, Phe290, and Phe331) stabilizes the methyl group of the bound ACh.[5] The catalytic mechanism of ACh hydrolysis is a successive two-step process, namely, acylation and deacylation steps, respectively. During hydrolysis, the imidazole ring of catalytic His440 acts as a general base to first increase the nucleophilicity of Ser200 hydroxyl oxygen in the acylation step and subsequently helps to activate a nucleophilic water molecule, which initiates the deacylation reaction.[6−15]

Figure 1

Illustration of (a) the overall AChE reactant structure and (b) the catalytic triad shown using ball-and-stick representation.

AChE is also reactive toward and is a primary target for an array of covalent inhibitors.[16−18] Specifically, the nucleophilic Ser200 of AChE is especially reactive toward organophosphate compounds; some examples of these compounds are illustrated in Figure S1 (Supporting Information). Organophosphate compounds prevent ACh hydrolysis by forming covalent adducts with the reactive Ser200 hydroxyl group, which results in overstimulation of the cholinergic synapses and leads to neuromuscular paralysis and fatality.[19−22] Nerve agents are military-grade organophosphorus compounds, whose synthesis began in the early 1930s and reached an all-time high during The Cold War.[23] Since then, nerve agents have been used in chemical warfare on multiple occasions and continue to pose a public health risk as potential chemical warfare agents.[23−26] Furthermore, less lethal organophosphate pesticides are widely used and unintentional exposure to these compounds represents a public health challenge.[21,22,27−30]

Soman (3,3-dimethylbutan-2-yl methylphosphonofluoridate) is a highly reactive nerve agent with a bimolecular inhibition rate constant of 0.86 ± 0.2 M–1 min–1 against huACHE.[31] Although computational studies of AChE phosphonylation[32−35] by other nerve agents such as sarin and tabun have been carried out, the AChE phosphonylation mechanism by soman has not been computationally characterized. Historically, the phosphonylation mechanism was assumed to proceed through an in-line displacement mechanism,[36−40] in which nucleophilic Ser200 attacks the phosphorus atom with concerted dissociation of the leaving group, as schematically shown in Figure 2a. However, recent studies suggested a stepwise addition–elimination mechanism,[13,32−35,41−46] in which the phosphonylation process begins with the formation of the phosphonate–ester covalent bond resulting in a trigonal bipyrimidal pentacovalent phosphonate intermediate, which collapses into the product by the dissociation of the leaving group; see Figure 2b. Alternately, phosphonylation could also proceed via a dissociative process, in which the fluorine dissociation precedes formation of the phosphonate–ester bond (Figure 2c). Several computational studies were carried out to study the AChE and nerve agent phosphonylation mechanism and revealed competing details regarding the rate-determining step,[32,34,35,46] the orientation of the leaving group with respect to the nucleophilic Ser200,[17,32,38,47,48] and last the role of catalytic Glu327 in phosphonylation.[33,35] Furthermore, no free energy profile has been determined for AChE phosphonylation reactions.

Figure 2

Three possible reaction pathways for AChE phosphonylation.

In the present study, we have characterized the phosphonylation reaction between AChE and soman and determined its free energy reaction profile for the first time. Our computational investigations center on Born–Oppenheimer ab initio QM/MM molecular dynamics simulations and the umbrella sampling method, a computational tour-deforce to study biochemical reactions:[14,49−61] the chemically reacting moieties are described by the ab initio QM method, the surrounding enzyme environment is treated explicitly by the classical MM force field, and the active site dynamics and those of the surroundings are simulated on an equal footing. This state-of-the-art computational approach has been demonstrated to be powerful in characterizing a number of complex systems, including the aging mechanism of soman inhibited AChE.[14,49−61]

Methods

The initial non-covalent AChE–soman complex was prepared on the basis of the crystal structure of non-aged TcAChE–soman (pdb code: 2WFZ).[4] On the basis of the crystal structure and biochemical results, the chiral phosphorus atom of soman was modeled to be in its S enantiomer. For each residue in the PDB file with partial occupation numbers, its configuration with the highest occupancy was retained. The protonation states of ionizable residues were determined by using H++[62] and by careful examination of the local hydrogen bond network. All histidine residues were modeled to be neutral with residues 440 and 486 protonated as HID and residues 26, 159, 181, 209, 264, 362, 398, 406, 425, and 471 protonated as HIE. Glu199 was modeled to be protonated as in our previous study,[61] while all other ionizable residues were modeled to be in their natural protonation states. Partial charges for soman were fitted with HF/6-31G(d) calculations and the restrained electrostatic potential (RESP) module in the Amber package and are listed in Table S1 (Supporting Information).[63] Other force field parameters were adapted from the GAFF[64] force field using Antechamber.[65]

AChE–P(S)-soman reactant was solvated to a rectangular box with a protein–box edge buffer distance of 10 Å and neutralized by adding counterions. The resultant model consisted of about 97 000 atoms and was simulated using the AMBER11[63] molecular dynamics package with the Amber99SB force field[66−68] and TIP3P[69] water model. The computational system was first energy minimized followed by a series of equilibration and a final 5 ns molecular dynamics production run with periodic boundary conditions. During the optimization and equilibration steps, the distance between catalytic Ser200 hydroxyl oxygen and soman phosphorus atoms was restrained to 3.2 Å and all restraints were released for the production run. Isotropic position scaling and the Berendsen thermostat method[70] were employed to maintain the system at 1 atm and 300 K, respectively. Long-range electrostatic interactions were treated with the particle mash Ewald (PME) method, and a 12 Å cutoff was used for both PME and van der Waals (vdW) interaction.

For the QM/MM model, a snapshot of the P(S)-soman in complex with AChE was taken from the 5 ns production molecular dynamics simulation. Ions and waters beyond 30 Å from the reaction center, chosen as the hydroxyl oxygen of Ser200, were deleted. The resultant QM/MM models had about 12 000 atoms, and the QM subsystem consisted of Ser200, His440, Glu327, and Tyr121 side chains and P(S)-soman, as shown in Figure S2 (Supporting Information). The QM/MM boundary was described by the pseudobond approach[50,71−75] with spherical boundary conditions. Only atoms within 23 Å of the reaction site were allowed to move, and 18 and 12 Å cutoffs were used for electrostatics and vdW interactions, respectively. No cutoff was employed for the electrostatic interactions between the QM and MM regions. The Beeman algorithm[76] was used to integrate Newton equations of motions using a time step of 1 fs. The Berendsen thermostat method[70] was employed to control the system temperature at 300 K. QM/MM simulations were run using modified versions of the Q-Chem[77] and Tinker[78] packages. The QM subsystem was described using the B3LYP functional with the 6-31G* basis set, and the MM subsystem was described using the Amber99SB force field and the TIP3P water model.

The QM/MM reactant models were first geometry optimized using the B3LYP/6-31G* level of theory, and several reaction paths have been examined using the reaction coordinate driving method.[50] In an effort to study whether phosphonylation is a concerted or stepwise reaction, two-dimensional minimal energy surfaces were mapped out along the addition and elimination reaction coordinates. For any given reaction path, the minimum energy path was first determined. Then, for every other configuration along the minimum energy path, a 500 ps MD simulation with a MM force field was performed to equilibrate the MM subsystem, while the QM subsystem was kept frozen. The resulting snapshots were used as the starting structures for Born–Oppenheimer QM/MM MD simulations with an ab initio QM/MM potential[50,79−83] and the umbrella sampling method.[84−86] Each QM/MM MD (umbrella window) simulation was run for 25 ps, and data along the last 20 ps were collected for analysis. The probability distributions along the reaction coordinates were then determined for each window and pieced together using the weighted histogram analysis method (WHAM).[87−89] The statistical error is estimated by averaging the free energy difference between 6 and 15 ps and 16 and 25 ps.

Results and Discussion

Figure 3a illustrates non-covalent interactions between P(S)-soman and the AChE active site: soman’s alkyl tail mainly interacts with aromatic residues of the AChE catalytic gorge, while soman’s phosphonyl oxygen forms hydrogen-bonding interactions with −NH groups of the enzyme’s oxyanion hole. In addition, Figure 3b schematically illustrates key residues of the active site that position soman in a reactive geometry within the active site. Hence, soman is positioned such that the nucleophilic oxygen of Ser200 can directly attack its phosphorus atom and initiate the irreversible inhibition process.

Figure 3

(a) Illustration of the AChE active site in the non-covalent complex with P(S)-soman, where the catalytic triad, oxyanion hole peptidic NH groups, and soman are shown using ball-and-stick representation. The frame corresponds to a randomly chosen snapshot collected from atomistic MD simulations and used as the starting reactant structure for QM/MM MD simulations. (b) Schematic illustration of structural interactions among AChE catalytic gorge residues and inhibitor soman; hashed lines represent hydrogen-bonding interactions, and the solid orange line represents the reactive path between AChE and soman.

The phosphonylation of AChE by soman involves two chemical events: one is the formation of a phosphonate–ester covalent bond between organophosphate and the hydroxyl oxygen of Ser200, and the other is the dissociation of the fluorine atom from the nerve agent. These are also termed addition and elimination steps, respectively. In order to examine their sequential order, first we mapped out a 2-D QM/MM minimal energy surface, as shown in Figure S3 (Supporting Information). The resultant surface clearly indicates that the phosphonate–ester bond formation precedes the dissociation of the fluorine atom, suggestive of an associative nucleophilic substitution mechanism. Snapshots along the identified minimum energy path, highlighted using dashed black lines in Figure S3 (Supporting Information), were used as the starting configurations to carry out MM equilibrations and ab initio QM/MM molecular dynamics simulations with umbrella sampling for the determination of the reaction free energy profile. As shown in Figure 4, this reaction free energy profile has two high-energy states that are separated by a shallow minimum corresponding to the trigonal bipyrimidal pentacovalent intermediate. The overall activation energy is ∼9.6 kcal/mol, and the rate-determining step corresponds to the dissociation of the fluorine atom. Configurations corresponding to stationary points along the phosphonylation reaction are illustrated in Figure 5, and the calculated key distances are shown in Figure S4 (Supporting Information).

Figure 4

Free energy profile for the addition–elimination mechanism of AChE phosphonylation by P(S)-soman determined along the last 20 ps of 25 ps QM/MM MD simulations. The reaction coordinate is defined as the linear combination of three bond lengths: RC = d(somanP···somanF) – d(Ser200HG···His440NE2) – d(Ser200OG···somanP). The total length of ab initio QM/MM MD simulations for this reaction path is 850 ps (34 windows × 25 ps each). The statistical error is estimated by averaging the free energy difference between 6 and 15 ps and 16 and 25 ps.

Figure 5

Illustration of key structures along the 63.1% associative phosphonylation reaction. The snapshots shown correspond to stationary points along the reaction, and the labeled distances (Å) were determined along the last 20 ps of the 25 ps QM/MM MD simulations.

In the reactant non-covalent complex (Figure 5 and Figure S5, Supporting Information), nucleophilic Ser200 is hydrogen bonded with the imidazole ring of catalytic His440 and this hydrogen bond helps increase the nucleophilicity of the enzyme. In this configuration, Ser200 hydroxyl oxygen is in an optimal orientation to attach the phosphorus atom of soman with an average P–OG distance of 3.21 ± 0.10 Å. The fluorine atom is axial with respect to the attacking nucleophile with an average OG–P–F angle of 166.02 ± 4.88°. Also, two of the three hydrogen bonds between the phosphonyl oxygen O1 of soman and the oxyanion hole are stable during QM/MM MD simulations. At the transition state (TS1) of the addition step, a partial phosphonate–ester bond forms between soman and AChE (average P–O distance of 2.10 ± 0.06 Å) and the Ser200 hydroxyl HG proton is shared between nucleophilic OG oxygen and NE2 nitrogen of the imidazole ring.

At the intermediate state (Figure 5 and Figure S6, Supporting Information), a metastable trigonal bipyrimidal pentacovalent phosphonate is formed and the third hydrogen bond between the oxyanion hole and phosphonyl oxygen also forms at this state (o1···o3). Here, the proton from the Ser200 hydroxyl group is completely transferred to His440 (average HG–NE2 distance of 1.05 ± 0.03 Å) and is hydrogen bonded to both OG and O2 oxygen atoms with an average H–O distance of 1.68 ± 0.13 and 2.44 ± 0.21 Å, respectively. The phosphorus atom is covalently bonded with five atoms (d1···d6), in which the phosphonyl oxygen (d4) maintains its double bonding characteristic, while the P–O (d1) and P–F (d3) bonds are longer than the P–O2 alkoxyl bond. In comparing the intermediate state with the reactant, we find the average P–F distance is increased by about 0.1 Å. The newly formed P–OG covalent bond is about 0.2 Å longer in the intermediate state compared with the product configuration. On the basis of Pauling’s formulation,[90,91] the characterized enzyme catalyzed phosphonylation reaction is 63.1% associative with an intermediate P–O distance of 1.85 ± 0.03 Å. In order to access the stability of the pentacovalent phosphonate, 28 snapshots were randomly selected around the intermediate state and used as the starting configurations for unrestrained QM/MM MD simulations. Within 1.5 ps, 92.9% of the simulations (26 out of 28) relaxed to the reactant state, while only 7.1% (2 out of 28) remained in the intermediate state configuration.

In the second elimination step, cleavage of the phosphorus–fluorine covalent bond results in a stable soman–AChE covalent conjugate and a negatively charged fluorine (F–) atom (Figure 5 and Figure S7, Supporting Information). Here, the hydrogen bonding interactions between the oxyanion hole and the phosphonyl O1 oxygen are maintained and help to stabilize the product. Meanwhile, the His440 imidazole ring is doubly protonated and its NE2 nitrogen is hydrogen bonded with OG and O2 oxygen atoms of the soman modified serine residue. The eliminated fluorine atom forms stabilizing interactions with Tyr121 hydroxyl hydrogen (f1), hydrogen atoms of soman’s alkyl chain (f3··· f6), and a water molecule (f2), which migrates into the functional site (Figure S7b, Supporting Information).

From the reactant to the product states, the mobility of the His440 imidazole ring is important for phosphonylation: His440 NE2 nitrogen is hydrogen bonded with the Ser200 hydroxyl group in the reactant state (average NE2–HG–OG angle of 164.30 ± 6.82°) and reorients during the reaction such that the NE2 nitrogen can form hydrogen bonds with soman O2 oxygen, connecting the phosphonate moiety to the akyl chain, in the product state (average NE2–HG–O2 angle of 156.40 ± 11.41°).

In order to access the validity of the characterized transition state (TS2), 36 snapshots were selected randomly around TS2 and used the starting structures for unrestrained QM/MM MD simulations. For each snapshot, two simulations were started with the Boltzmann sampled and reversed velocities. Out of the 36 pairs of simulations, 6 pairs relaxed toward the reactant state and 13 pairs relaxed toward the product, while 17 pairs relaxed either forward toward the product or backward toward the reactant. Overall, 8.3% of the simulations converged to the reactant state, while 51.4% converged to the intermediate and 40.3% of the simulations converged to the product state, suggesting that the characterized transition state (TS2) is meaningful. As shown in Figure 6, we can see that the active site of the computationally determined non-aged soman–AChE conjugate is very similar to the experimentally determined crystal structure configuration, which provides further support to our characterized mechanism.

Figure 6

Superimposition of the active site of non-aged soman-inhibited AChE crystal structure[4] and the phosphonylated product determined from QM/MM MD simulations, shown using cyan and green carbons, respectively.

Conclusions

By employing Born–Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, a state of the art computational approach to study chemical reactions in biological systems, we have elucidated the irreversible inhibition/phosphonylation reaction mechanism of AChE and P(S)-soman and determined its free energy profile for the first time. It should be noted that there is no previous theoretical study on AChE phosphonylation by soman, and its mechanism has been ambiguous so far, with three possible reaction schemes suggested in the literature. Our computational results revealed a highly coupled associative addition–elimination reaction mechanism for this covalent inhibition process. During the enzyme-catalyzed reaction, the addition and elimination steps are separated by a metastable pentacovalent phosphonate intermediate. Catalytic His440 of AChE initiates the phosphonylation reaction by deprotonating the nucleophilic Ser200. Ser200 side chain oxygen can then attack soman’s phosphorus. The leaving fluorine atom is located in an axial orientation with respect to the attacking group. Together with our previous characterization of the aging mechanism of soman inhibited AChE, our simulations have revealed new detailed mechanistic insights into the damaging function of the nerve agent soman.