You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


1. Introduction {#sec1}
===============

The physiological role of acetylcholinestrase (AChE) is to effectively terminate cholinergic signals by hydrolyzing acetylcholine (ACh).^[@ref1],[@ref2]^ AChE has a general α/β hydrolase globular fold, as illustrated in Figure [1](#fig1){ref-type="fig"}a. 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 [1](#fig1){ref-type="fig"}b).^[@ref3]^ (*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,^[@ref3],[@ref4],[@ref23],[@ref26]^ while the acyl pocket (Trp233, Phe228, Phe290, and Phe331) stabilizes the methyl group of the bound ACh.^[@ref5]^ 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.^[@ref6]−[@ref15]^

![Illustration of (a) the overall AChE reactant structure and (b) the catalytic triad shown using ball-and-stick representation.](jp-2014-02712d_0001){#fig1}

AChE is also reactive toward and is a primary target for an array of covalent inhibitors.^[@ref16]−[@ref18]^ Specifically, the nucleophilic Ser200 of AChE is especially reactive toward organophosphate compounds; some examples of these compounds are illustrated in Figure S1 ([Supporting Information](#notes-1){ref-type="notes"}). 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.^[@ref19]−[@ref22]^ Nerve agents are military-grade organophosphorus compounds, whose synthesis began in the early 1930s and reached an all-time high during The Cold War.^[@ref23]^ 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.^[@ref23]−[@ref26]^ Furthermore, less lethal organophosphate pesticides are widely used and unintentional exposure to these compounds represents a public health challenge.^[@ref21],[@ref22],[@ref27]−[@ref30]^

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 *hu*ACHE.^[@ref31]^ Although computational studies of AChE phosphonylation^[@ref32]−[@ref35]^ 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,^[@ref36]−[@ref40]^ in which nucleophilic Ser200 attacks the phosphorus atom with concerted dissociation of the leaving group, as schematically shown in Figure [2](#fig2){ref-type="fig"}a. However, recent studies suggested a stepwise addition--elimination mechanism,^[@ref13],[@ref32]−[@ref35],[@ref41]−[@ref46]^ 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 [2](#fig2){ref-type="fig"}b. Alternately, phosphonylation could also proceed via a dissociative process, in which the fluorine dissociation precedes formation of the phosphonate--ester bond (Figure [2](#fig2){ref-type="fig"}c). Several computational studies were carried out to study the AChE and nerve agent phosphonylation mechanism and revealed competing details regarding the rate-determining step,^[@ref32],[@ref34],[@ref35],[@ref46]^ the orientation of the leaving group with respect to the nucleophilic Ser200,^[@ref17],[@ref32],[@ref38],[@ref47],[@ref48]^ and last the role of catalytic Glu327 in phosphonylation.^[@ref33],[@ref35]^ Furthermore, no free energy profile has been determined for AChE phosphonylation reactions.

![Three possible reaction pathways for AChE phosphonylation.](jp-2014-02712d_0002){#fig2}

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:^[@ref14],[@ref49]−[@ref61]^ 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.^[@ref14],[@ref49]−[@ref61]^

2. Methods {#sec2}
==========

The initial non-covalent AChE--soman complex was prepared on the basis of the crystal structure of non-aged *Tc*AChE--soman (pdb code: 2WFZ).^[@ref4]^ 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++^[@ref62]^ 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,^[@ref61]^ 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](#notes-1){ref-type="notes"}).^[@ref63]^ Other force field parameters were adapted from the GAFF^[@ref64]^ force field using Antechamber.^[@ref65]^

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^[@ref63]^ molecular dynamics package with the Amber99SB force field^[@ref66]−[@ref68]^ and TIP3P^[@ref69]^ 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^[@ref70]^ 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](#notes-1){ref-type="notes"}). The QM/MM boundary was described by the pseudobond approach^[@ref50],[@ref71]−[@ref75]^ 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^[@ref76]^ was used to integrate Newton equations of motions using a time step of 1 fs. The Berendsen thermostat method^[@ref70]^ was employed to control the system temperature at 300 K. QM/MM simulations were run using modified versions of the Q-Chem^[@ref77]^ and Tinker^[@ref78]^ 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.^[@ref50]^ 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^[@ref50],[@ref79]−[@ref83]^ and the umbrella sampling method.^[@ref84]−[@ref86]^ 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).^[@ref87]−[@ref89]^ The statistical error is estimated by averaging the free energy difference between 6 and 15 ps and 16 and 25 ps.

3. Results and Discussion {#sec3}
=========================

Figure [3](#fig3){ref-type="fig"}a 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 [3](#fig3){ref-type="fig"}b 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.

![(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.](jp-2014-02712d_0003){#fig3}

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](#notes-1){ref-type="notes"}). 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](#notes-1){ref-type="notes"}), 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](#fig4){ref-type="fig"}, 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](#fig5){ref-type="fig"}, and the calculated key distances are shown in Figure S4 ([Supporting Information](#notes-1){ref-type="notes"}).

![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*(soman~P~···soman~F~) -- *d*(Ser200~HG~···His440~NE2~) -- *d*(Ser200~OG~···soman~P~). 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.](jp-2014-02712d_0004){#fig4}

![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.](jp-2014-02712d_0005){#fig5}

In the reactant non-covalent complex (Figure [5](#fig5){ref-type="fig"} and Figure S5, [Supporting Information](#notes-1){ref-type="notes"}), 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](#fig5){ref-type="fig"} and Figure S6, [Supporting Information](#notes-1){ref-type="notes"}), 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,^[@ref90],[@ref91]^ 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](#fig5){ref-type="fig"} and Figure S7, [Supporting Information](#notes-1){ref-type="notes"}). 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](#notes-1){ref-type="notes"}).

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](#fig6){ref-type="fig"}, 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.

![Superimposition of the active site of non-aged soman-inhibited AChE crystal structure^[@ref4]^ and the phosphonylated product determined from QM/MM MD simulations, shown using cyan and green carbons, respectively.](jp-2014-02712d_0006){#fig6}

4. Conclusions {#sec4}
==============

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.

Figures showing the chemical structures of some organophosphate compounds; the QM/MM partition for the non-covalent AChE--soman reactant complex; the two-dimensional minimal energy surface of the AChE--soman phosphonylation reaction; calculated key distances during the characterized addition--elimination phosphonylation reaction mechanism; illustrations of the non-covalent reactant complex, the trigonal bipyrimidal pentacovalent phosphorus intermediate, and the covalent AChE--soman conjugate product; and distances in QM/MM MD simulations. Tables showing partial charges, atom types, and connectivity assignments for soman and some of the key geometric parameters for the phosphonylation/inhibition reaction of AChE--soman. This material is available free of charge via the Internet at <http://pubs.acs.org>.

Supplementary Material
======================

###### 

jp502712d_si_001.pdf

The authors declare no competing financial interest.

This work was supported by NSF-IGERT and NIH (R01-GM079223). We thank NYU-ITS and NYUAD for providing computational resources.