Epoxide hydrolysis as a model system for understanding flux through a branched reaction scheme

Kinetic rates and computer simulations can explain regioselectivity and stereoconfiguration in the products of enzyme-catalyzed epoxide hydrolysis.


S1. Simulation Details
Following from our previous work on the StEH1-catalyzed hydrolysis of styrene oxide and trans-stilbene oxide, 1,2 we have used the empirical valence bond (EVB) approach 3 to model the chemical step for the hydrolysis of 1a. The substrate and key intermediates were parameterized as in our previous work, using the OPLS-AA force field, 4

version 11, as implemented in
Macromodel v. 9.1. 5 Partial charges were obtained using the restrained electrostatic potential (RESP) approach, 6 by calculating the electrostatic potential at the HF/6-31G* level of theory, using Gaussian 09 Rev. C.01, 7 and fitting the potential to the atomic coordinates using ANTECHAMBER. 8 All EVB parameters used in this work are presented in as Supplementary Tables below. All molecular dynamics (MD) and EVB simulations were performed using the OPLS-AA force field 3 and the Q simulation package. 9 All energetic analysis was performed using the Qfep and Qcalc modules of Q, and all geometric analysis was performed by postprocessing of the Q trajectories using GROMACS v. 5.0.5. 10 Coordinates for the wild-type enzyme and the C1 variant were taken from the Protein Data Bank, 11 PDB IDs 2cjp 12 and 4uhb 13 respectively. We note that the K141 side chain (V141K substitution from wild-type) is poorly resolved in the C1 structure; however, as this is a solvent exposed lysine that is likely to be highly flexible, we kept the crystallographic coordinates in our simulations. As in our previous work, 1 both enantiomers of 1a were manually placed in the active site in such a way as to optimize the angle of attack for the nucleophile, protein-substrate interactions, and prevent unbinding of the substrate during the simulations. Each enantiomer was placed in each of the two binding modes shown in Figure 2 of the main text, leading to a total of four independent enantiomer/binding mode combinations per enzyme variant.
The structures were prepared for the EVB simulations as in our previous studies. 1,2 In brief, the system was solvated by a 20Å droplet of TIP3P water molecules, 14 centered on the γ-carbon of S3 D105, and with the bonds to hydrogen atoms restrained using the SHAKE algorithm. 15 This droplet was then described using a multi-layer model, in the framework of the surface constrained all atom solvent approach (SCAAS) as implemented in Q v. 5.0, 9 in which all atoms in the inner 85% of the sphere (17Å from the sphere center) were unrestrained and fully flexible during our simulations, all atoms in the outer 15% of the sphere (outer 3Å layer) were subjected to 20 kcal mol -1 Å -2 harmonic restraint to reduce their mobility, and all atoms outside the explicit sphere were restrained to their crystallographic positions using a 200 kcal mol -1 Å -2 harmonic restraint. Tying in with this, all amino acid side chains within the mobile region, i.e. the inner 17Å of the sphere, were ionized according to the most likely ionization state predicted by PROPKA (v. 3.0), 16,17 and all residues located outside this sphere were kept in their neutral form to avoid unsolvated charges from exerting unphysical electrostatic effects on the reacting region, in line with our previous work. 1,18,19 The most likely histidine protonation states were predicted using the MOLPROBITY server, 20 and H104 was kept protonated, in line with our previous work. 1,2 This led to an overall system charge of -1 for both variants. For a list of residues ionized during the simulations, and histidine protonation states, see the supporting information of ref. 2.
Once the structures had been prepared for the simulations, each system (at the Michaelis complex) was first minimized at close to 0K with an extremely small step size to ease initial bad contacts, after which the step size was set to 1 fs for the remainder of the simulation, with a final MD minimization at 1K for 30 ps. The system was equilibrated by gradually heating the water sphere for 30 ps at 100K and another 30 ps at 300K, while keeping all solute atoms (including the catalytic water and substrate) restrained to their crystallographic positions with a 200 kcal mol -1 Å -2 harmonic restraint. The system was then cooled down again to 5K, and then gradually reheated to first 30, then 150, then 300K, while dropping the restraints from 200 kcal mol -1 Å -2 on all solute atoms to 200 kcal mol -1 Å -2 on only the reacting atoms ( Figure S1). Once this heating-S4 cooling-heating procedure was complete, the system was then equilibrated for a further 10 ns as preparation for our EVB simulations, and the root mean square deviations (RMSD) of all C α atoms during out simulations are shown in Figure S3 for the wild-type and R-C1 variant respectively, and the corresponding values for the substrate are shown in Figure S4. A 1 fs step size was used throughout, and temperature was kept constant using Berendsen's thermostat 21 with a bath coupling of 100 fs. During the whole procedure, the non-bonded pair lists were updated every 30 ps, with a cutoff for the explicit interactions at 10 Å, while long range electrostatic interactions beyond this distance were modelled using the Local Reaction Field (LRF) 22 approach without any cutoff. The interactions of the atoms in the reacting region had no cutoff assigned to them and explicitly interacted with the complete system. This procedure was repeated 30 times per system, using different initial velocities, to a total simulation time of 300 ns per system and 2.4 μs over all systems. In each case, the final equilibrated structure was then used as the starting point for 30 independent EVB simulations per system. As in our previous work, all EVB calculations were performed in a two-step process, with the reacting system described using the valence bond states shown in Figure S1. The EVB free energy perturbation/umbrella sampling (FEP/US) calculations were performed in 51 individual mapping windows of 200 ps in length each, to a total simulation time of 10.2 ns per system and 306 ns per reacting system/trajectory (4.9 μs over all trajectories). The EVB simulations were performed in two sequential steps, simulating first nucleophilic attack on the epoxide ring to yield the alkylenzyme intermediate, followed by water attack on the alkylenzyme intermediate to yield the transition state for the formation of the tetrahedral intermediate, using the end-points of the first simulation step as starting points for the second simulation step (see the main text for the reaction mechanism). From this point on, the velocities of all atoms in the systems were regenerated according to the Maxwell distribution at 300 K, before letting the system follow the S5 path towards both the alkylenzyme and tetrahedral intermediates again, using the same amount of sampling as for the first reaction step. This was necessary to obtain the correct alignment of the nucleophilic water molecule for proton abstraction by H300 during the formation of the tetrahedral intermediate. System energies were saved every 10 fs, and used to generate the EVB free energy profiles.
The corresponding off-diagonal and gas-phase shift parameters necessary to generate the EVB free energy profiles for each system are presented in Section S5 (for details of what these parameters are, see for example reviews in refs. 23 and 24. These parameters were obtained by modeling the corresponding uncatalyzed reaction in aqueous solution using the substrate, nucleophilic water molecule, methyl imidazole and propionate (as models for the D105/H300 side chains, respectively), and fitting the free energy profiles for the corresponding uncatalyzed reaction to results obtained from quantum chemical calculations. The quantum chemical calculations were performed using Gaussian 09, revision E.01, with the system described using the ωB97Χ-D functional. Stationary points were optimized using the 6-31+G* basis set and impicit solvation using the SMD method, while frequency calculations were performed at the 6-311+G** level of theory. The intrinsic reaction coordinates were followed from the transition state to the respective reactant and product states, with the final energies calculated for the fragments at infinite separation. All pathways were calculated individually, with the lowest energy one then being used to parametrize the respective reference reactions.
For the quantum calculations, the transition states obtained from the EVB calculations were further optimized using the functional mentioned above towards the stationary points. From this point, the internal reaction coordinates (IRC) 25,26 were followed both towards the reactant and product states to obtain the QM potential energy surface for the first reaction step of the enzymatic mechanism. The final points of the IRC were again optimized to the stationary points, S6 together with the reactants at infinite separation. The resulting energetics were used for fitting the EVB parameters in aqueous solution for each system (attack at C-1 and C-2 for each enantiomer; note that the two binding modes are expected to have the same energy in aqueous solution as the environment is homogenous). The corresponding results are presented in Table   S1. All parameters were then transferred unchanged to the corresponding enzymatic system, which is made feasible by the phase-independence of the EVB off-diagonal term as discussed in detail in for example ref. 27.
Finally, Figures S5 and S6 show the structures of key stationary points for the hydrolysis of (R,R)-and (S,S)-1a catalyzed by wild-type StEH1 respectively, and the corresponding geometric parameters are shown in Table S2.

S2. Supplementary Scheme
Scheme S1. Model of kinetic mechanism for epoxide hydrolysis catalyzed by StEH1. The model is based on data from analysis of pre-steady state and steady state kinetics as well as the enantiomeric ratios of the product diols.    S12 Figure S5: Key stationary points for the hydrolysis of (R,R)-1a by wild-type StEH1. Coordinates correspond to clustered structures obtained from trajectories at the respective extrema of the EVB free energy profiles. S13 Figure S6: Key stationary points for the hydrolysis of (S,S)-1a by wild-type StEH1. S14 Table S1: A comparison of the energetics of the uncatalyzed hydrolysis of both (R,R)-and (S,S)-MeSO, from our DFT calculations and EVB calibration procedure a . The corresponding absolute energies for the DFT calculations are shown in Table S2. 10.0 ± 0.03 10.0 ± 0.02 10.0 ± 0.03 10.0 ± 0.01 a DFT calculations performed as described in Section S2. In the second step, which involves ester hydrolysis by nucleophilic attack of a water molecule on an oxirane intermediate (Figure S2), all EVB energetics were calibrated to a ∆G ‡ of 19 kcal mol -1 and a ∆G0 of 10 kcal mol -1 , as described in Section S1 and the Supporting Information of our previous work. 1 As the two binding modes described in     Step  a For all atoms except reacting atoms, a standard 6-12 Lennard Jones potential was used. In the case of the reacting atoms, which change bonding patterns between atoms i and j, an alternate function of the form: Vreact = Ci Cj exp (-αiαjrij) was used to prevent artificial repulsion between these atoms as bonding patterns change. rij denotes the distance (Å) between atoms i and j.         a Torsion angle potential: Vtorsion= V1 (1+cos(nφ -δ))+ V2 (1+cos2(nφ -δ))+ V3 (1+cos3(nφ -δ)), n is the periodicity (number of maxima per turn) and δ is the phase shift.

Improper Type
Not Set S34 Table S15: Improper torsion angles of the VB states used to describe the reacting system, following the atom numbering shown in Figure S7.  were microsolvated with two water molecules, and each of the reacting complexes were microsolvated with four water molecules.