Biochemistry other enzymes.4,39,56−58 Further details of the setup of each model system, including addition of hydrogen atoms, solvation, minimization, and stochastic boundary conditions, are provided in the Supporting Information. Initial QM/MM MD Simulations. Starting structures for the QM/MM umbrella sampling MD simulations were generated by performing ten 14 ps stochastic boundary QM(AM1)/MM(CHARMM22) MD simulations on the same starting structure, for each enzyme−substrate complex. Each simulation was given a different set of initial velocities by supplying a different seed value to the random number generator, to generate different starting conformations for the umbrella sampling MD simulations. For each enzyme−substrate complex, the substrate and nucleophilic Asp333 side chain were treated quantum mechanically at the AM1 level of theory.59 The rest of the enzyme and solvent were treated with the CHARMM22 MM force field.60 A hydrogen link atom61 was used to partition the covalent bond between the QM and MM regions, placed between the Cβ and Cα atoms of the Asp333 residue. QM/MM Umbrella Sampling MD Simulations. Ten simulations of the nucleophilic attack at C(1) or C(2) of the epoxide ring were conducted for each substrate (defined in Figure 1) using QM(AM1)/MM(CHARMM22) umbrella sampling MD. The coordinates of the final structures from the unrestrained simulations described above were used as the initial structures. Multiple starting geometries were used to enhance the amount of conformational sampling, which has been found to be important in other enzyme studies.4 On the basis of each starting geometry, the reaction was modeled using umbrella sampling with conditions and protocols similar to those developed and tested for QM/MM simulations of other enzymes.40,56,57,62 The reaction coordinate (denoted herein as r) was defined as the difference between two interatomic differences: the epoxide O−C(1) or −C(2) distance minus the Oδ2 (Asp333)−C(1) or −C(2) distance, corresponding to the bonds being broken or formed, respectively. Similar reaction coordinates have been used in previous studies of other enzymes.40,56,57 The Oδ1 atom of Asp333 is not in a position to favor nucleophilic attack at either epoxide carbon atom; hence, only nucleophilic attack by Oδ2 is considered. Series of simulations were performed with the reaction coordinate harmonically restrained at values between −2.0 and 0.5 Å in 0.05 Å steps, using a force constant of 100 kcal mol−1 Å−2. All other aspects of the umbrella sampling MD simulations were identical to those used for the unrestrained MD simulations, described above. Each simulation consisted of 5 ps of equilibra- tion dynamics and 25 ps of sampling dynamics, initiated by using the end point of the first 5 ps from the previous simulation.58,63 The weighted histogram analysis method (WHAM)64 was used to combine the reaction coordinate statistics of the simulations along the reaction coordinate. Multiple umbrella sampling MD runs were performed, to test the validity and reproducibility of the results. The free energy barriers were calculated by combining the statistics from all of the simulations performed for each reaction in the WHAM unbiasing calculation. Such approaches have been applied successfully in previous work for other enzymes.40,56,57 Preliminary umbrella sampling simulations were performed with the His523 Nε atom either protonated or unprotonated, to identify the most likely protonation state of this residue. In these simulations, five different simulations were run using the same setup and QM/MM methodology described above, but with the following exceptions: r was varied in steps of 0.1 Å, Article and the simulations consisted of 1 ps of equilibration dynamics and 9 ps of sampling dynamics. Higher-Level QM/MM Calculations and Energy Cor- rections. AM1, like other semiempirical QM methods, often overestimates activation barriers for chemical reactions (by ≥10 kcal/mol).53,65 Higher-level QM methods provide more accurate calculations of activation barriers, but their computa- tional cost prohibits MD sampling. Instead, the AM1-based free energy barriers were corrected afterward using the following (previously used) protocol.65 The end point of the umbrella sampling simulation with the value of r corresponding to the highest point on each free energy surface (i.e., corresponding to the transition state) was used as the initial geometry for AM1/ CHARMM22 potential energy calculations. Potential energy surfaces (PESs) were calculated in CHARMM: a series of geo- metry optimizations was performed, with the reaction coordinate (r) fixed successively to a range of values. These PESs were initiated at an r of −0.2 Å (transition state) and followed downhill in energy to r values of −2.0 and 0.5 Å, corresponding to the reactant complex and covalent ester intermediate, respectively (details in the Supporting Information). As with all other energy barriers calculated in this work, energy barriers were calculated as the difference in energy between the transition state and reactant complex. An activation entropy65 was calculated by computing the difference between the free energy barrier mentioned above and the Boltzmann-weighted average of the AM1/CHARMM22 potential energy barriers (ΔEa⧧ve) for the same process, using the expression in eq 1.4 (1) ⧧ ΔEave = −RT ln ∑ n i=1 ⎛ ⧧ ⎞ exp⎜⎜ −ΔEi ⎟ ⎝ RT ⎟⎠ The activation entropy was added to the B3LYP/CHARMM22 and SCS-MP2//B3LYP/CHARMM22 barriers (detailed below) to yield approximate free energy barriers. This approach gives only an approximate estimate of the activation entropy but has been shown to be reasonable for other enzymes. Composite calculations of this sort have previously given activation free energies in “quantitative” agreement with experiment for other enzymes, such as chorismate mutase and p-hydroxy- benzoate hydroxylase.65,66 The activation entropy is small for sEH and is not expected to differ significantly between different substrates; hence, the same average entropy correction is applied in all models. Entropic contributions to activation free energy barriers have been calculated for the enzyme subtilisin in a more rigorous manner by Warshel et al.67 and were shown to be small (∼2.5 kcal/mol at 300 K). PESs were calculated with (DFT)QM/MM methods using QoMMMa,68 which uses Tinker69 with the CHARMM2260 force field for the MM part of the calculation and Jaguar70 for the QM part. PESs for all systems were calculated at the B3LYP/6- 31G(d)/CHARMM22 level using the transition state structures from the corresponding AM1/CHARMM22 PESs as initial geo- metries. The reaction coordinate used in the AM1/CHARMM22 calculations was also used for the DFT/MM calculations. Structures were fully optimized at the DFT level, subject to reaction coordinate constraints. All atoms located more than 20 Å from the center of the system (Oδ2 of Asp333) were held fixed. Two PESs for t-SO were also optimized at the BP86/6-31G(d)/ CHARMM22 and BHandHLYP/6-31G(d)/CHARMM22 levels of theory, to test the dependence of the barriers on the den- sity functional. The basis set dependence of the B3LYP/ CHARMM22 energies was investigated by performing 1777 dx.doi.org/10.1021/bi201722j | Biochemistry 2012, 51, 1774−1786