Understanding the Molecular Mechanism of the Ala-versus-Gly Concept Controlling the Product Speci ﬁ city in Reactions Catalyzed by Lipoxygenases: A Combined Molecular Dynamics and QM/MM Study of Coral 8 R ‑ Lipoxygenase

: Lipoxygenases (LOXs) are a family of enzymes that catalyze the highly speci ﬁ c hydroperoxidation of polyunsaturated fatty acids, such as arachidonic acid. Di ﬀ erent stereo- or/and regioisomer hydroperoxidation products lead later to di ﬀ erent metabolites that exert opposite physiological e ﬀ ects in the animal body and play a central role in in ﬂ ammatory processes. The Gly-Ala switch of a single residue is crucial for the stereo- and regiocontrol in many lipoxygenases. Herein, we have combined molecular dynamics simulations with quantum mechanics/molecular mechanics calculations to study the hydrogen abstraction step and the molecular oxygen addition step of the hydroperoxidation reaction of arachidonic acid catalyzed by both wild-type Coral 8 R -LOX and its Gly427Ala mutant. We have obtained a detailed molecular understanding of this Ala-versus-Gly concept. In wild type, molecular oxygen adds to C 8 of arachidonic acid with an R stereochemistry. In the mutant, Ala427 pushes Leu385, blocks the region over C 8 , and opens an oxygen access channel now directed to C 12 , where molecular oxygen is added with an S stereochemistry. Thus, the speci ﬁ city turns out to be dramatically inverted. Since Leu385 is highly conserved among many lipoxygenase isoforms, this mechanism can be general, and we propose that the presence of such type of bulky and hydrophobic residues can be key in controlling the extreme regio- and stereospeci ﬁ city of lipoxygenases and, as a consequence, their physiological e ﬀ ects. Scheme of the reaction mechanism, arachidonic acid structure, MD simulations in WT Coral 8 R -LOX, hydrogen abstraction step in WT Coral 8 R -LOX, oxygen insertion step in WT Coral 8 R -LOX, conserved residues of the active site cavity of 8 R -LOX, free energy calculations of the oxygen insertion step in WT Coral 8 R -LOX, MD simulations in the mutant Gly427Ala 8 R LOX, oxygen insertion step in the mutant Gly427Ala 8 R LOX, free energy calculations of the oxygen insertion step in the mutant Gly427Ala 8 R -LOX, Leu385, and Leu431


INTRODUCTION
Most chemical catalysts are not very selective and catalyze a broad set of reactions. Conversely, enzymes are usually highly selective, catalyzing only specific reactions. This specificity is many times, due to the very subtle, small molecular details of the enzymes. A deep understanding of the way the enzymes govern the regio-and stereospecificity of the reactions they catalyze is key to design new enzymes with predetermined functions. Especially interesting is the case of lipoxygenases (LOXs). Lipoxygenases are a family of non-heme iron dioxygenases that catalyze the highly specific hydroperoxidation of polyunsaturated fatty acids containing one or several 1,4-Z,Zpentadiene units. 1−5 Animal LOXs are named depending on the specific position of oxygenation of their main substrate, arachidonic acid (AA; C20:4, ω-6), a carbon tetra-unsaturated fatty acid. 4,6−10 For instance, human 15-LOX, human 12-LOX, mouse 8-LOX, and human 5-LOX very specifically first abstract a hydrogen atom (the rate-determining step of the global process) from C 13 , C 10 , C 10 , or C 7 (a bisallylic methylene carbon) of AA, respectively (see Scheme S1). Subsequently, in a second step, molecular oxygen is predominantly added to C 15 , C 12 , C 8 , or C 5 of AA, respectively. This is followed by a retrohydrogen transfer to the formed peroxy radical to give the final hydroperoxy product, a hydro(pero)xyeicosatetraenoic acid (H(p)ETE). 3,11 This regioselectivity is accompanied in each case by the corresponding high stereoselectivity (S or R).
Key mutations have proven to radically modify the specificity of some lipoxygenases. For instance, whereas wild-type human 5-LOX almost exclusively oxygenates position 5 of arachidonic acid (with an S stereochemistry), the quadruple mutant F359W/A424I/N425M/A603I 5-LOX exhibits a major 15Slipoxygenase activity (85−95%). 9,12 What particular mutations have to be introduced to produce a concrete change of the enzyme functionality depends on each case. However, Coffa and Brash have identified a single residue that is crucial for stereo-and regiocontrol in many lipoxygenases: it is conserved as an Ala in S lipoxygenases, but as a smaller Gly in R lipoxygenases. 6,8,13,14 A Gly-Ala mutagenetic switch changes the oxygenation position in arachidonic acid, along with R or S stereochemistry (see Figure S1). 6 For instance, human 15-LOX-2 exclusively gives a 15S oxygenation product, whereas its Ala416Gly mutant mainly leads to 11R-H(p)ETE, with lower amounts of 15S-H(p)ETE. Conversely, coral 8R-LOX produces only 8R-H(p)ETE, but its Gly427Ala mutation almost exclusively gives 12S-H(p)ETE. Similar alterations have been observed in another R-LOX, human 12R-LOX, in which the Gly411Ala mutant produces the 8-HETE S-stereoisomer. All these results have allowed the formulation of the so-called Alaversus-Gly concept. Recent mutagenesis experiments 15 on a wide variety of S-lipoxygenases have shown that an Ala-to-Gly exchange does not necessarily lead to a complete inversion of the reaction specificity, although in most cases the formation of a significant amount of the R-oxygenation product is induced by that mutation.
On the other hand, the number of S-LOXs in nature is much greater than their counterpart R-LOXs, and much less is known about the latter ones. Coral (Plexaura homomalla) 8R-LOX is one of the family members whose crystallographic structure is solved. Its sequence is around 40% identical to that of human 5-LOX. 16 The crystal structure of the ligand-free form of 8R-LOX was reported a few years ago (PDB ID 3FG1), 17 and more recently, a complete structure of the 8R-LOX:AA complex was released by the same authors (PDB ID 4QWT). 18 In this structure, unlike the case of human 5-LOX, 19 the substrate coordinates are accurately determined, and a tail-first substrate orientation with the AA carboxylate end interacting with Arg182 ( Figure 1) is observed.
Experimentally it has been observed that the hydrogen abstraction step by coral 8R-LOX comes essentially from a unique position, pro-R H 10 . 18 In the crystallographic structure of the AA-bound form, the substrate is well-positioned to yield the hydrogen abstraction from C 10 ( Figure 1). More interesting though is the second reaction step, the oxygen insertion. As seen above, coral 8R-LOX presents a paradigm for the Alaversus-Gly concept.
The complete molecular basis of the Ala-versus-Gly concept is not well understood yet. In coral 8R-LOX Newcomer and coworkers 17 have suggested that the invariant Leu431, which is spatially close to Gly427 (one turn of the helix far), can be involved in controlling the oxygen access, but this assumption is based only on the analysis of the crystallographic structure. They propose that the wild-type enzyme (containing Gly427) leads exclusively to the 8R product because Leu431 shields C 12 from oxygenation. The Gly427/Ala427 switch, so introducing an additional methyl group, would induce a displacement of Leu431, pushing it away from C 12 , but now hindering the oxygen access to C 8 , thus producing the 12S product.
From the theoretical point of view, to the best of our knowledge, there is only one work describing the mechanism of coral 8R-LOX to date. 20 However, the authors modeled the substrate arachidonic acid as the quite shorter (4Z,7Z)-1,4,7,10undecatetraene, which was added in the active site manually, only carried out an extremely short molecular dynamics simulation of 100 ps, and did not study the O 2 addition to C 12 . Taking into account that an accurate description of the positioning of the actual substrate inside the cavity of the enzyme is absolutely crucial to account for the full catalytic mechanism of the enzyme (including its regio-and stereospecificity), and that this positioning highly depends on a complicated interplay of stabilizing interactions and steric hindrance between the actual substrate and the lipoxygenase, further work is needed to discuss the factors that rule the specificity of the catalytic reaction.
The purpose of this paper is 2-fold. First, we study theoretically the mechanism of the hydroperoxidation of arachidonic acid catalyzed by the enzyme coral 8R-LOX, one of the few reported R lipoxygenases so far. A number of papers have been devoted to the theoretical study of the catalytic mechanism of several S lipoxygenases by modeling the entire system, 21−32 but this is the first time that a similar study has been done for an R lipoxygenase. Second, we uncover insights on the factors that govern the regio-and stereospecificity of this reaction. Thus, a main goal will be to provide a detailed molecular understanding of the Ala-versus-Gly concept. To this aim, we have combined molecular dynamics (MD) simulations with quantum mechanics/molecular mechanics (QM/MM) calculations to study the hydrogen abstraction step and the molecular oxygen addition step of the hydroperoxidation reaction of arachidonic acid catalyzed by both wild-type coral 8R-LOX and its Gly427Ala mutant.

METHODOLOGY
2.1. Preparation of the Initial Coordinates of the 8R-LOX:AA Complex. Crystal structure 4QWT of 8R-LOX contains the AA substrate in subunit C. This subunit together with its substrate has been selected and prepared to study the reactivity of this system. The crystallographic coordinates lack a small region (indexes 307−315) in all four subunits. This loop region has been reconstructed in subunit C by the "loop refinement" facility of the program Modeler 9.2. 33 2.2. In Silico Mutagenesis. The starting structure employed to generate the in silico mutant Gly427Ala corresponds to the hydrogen abstraction product of the WT 8R-LOX (see Results), as this mutation has been proposed to affect only the second reaction step. The mutant was generated

ACS Catalysis
Research Article by changing the residue Gly427 to Ala with the Mutator Plugin of the VMD program. 34 2.3. Tunnel Representations. The oxygen access channels have been represented with the Caver3.0 program 35 used as a plugin of the PyMOL visualization program. In the WT system, the C8 carbon atom was used as initial starting point for the tunnel computation, whereas, in the Gly427Ala mutant, the C12 carbon atom was selected as the starting point. The default tunnel computation settings were selected, except for the minimum probe radius, for which a value of 0.75 Å was used.
2.4. Molecular Dynamics Simulations. All hydrogen coordinates have been generated with the program PROPKA version 3, 36,37 using a standard pH = 7 for the tritable residues. The system has been solvated in a cubic box of pre-equilibrated TIP3P water molecules of dimensions 125 Å × 125 Å × 125 Å. Then the system was neutralized with K + ions, and an additional KCl salt concentration was added to reproduce the experimental ionic strength of a 0.15 M solution. The resulting system contains nearly 184000 atoms, of which about 11000 belong to the protein.
MD simulations were run with the program NAMD2.9. 38 The CHARMM36 force field 39 was used for the protein and lipid atoms, and the Fe coordination sphere was described by an update of the parameters specifically derived by Saam et al. 40 First, the system was sent to 10000 initial energy minimization steps using the Conjugate Gradient algorithm. Then, MD simulations using periodic boundary conditions (PBC) were carried out. An equilibration step of 1 ns at constant temperature and pressure (303.15 K and 1 atm) was run. Langevin dynamics was used to control temperature, and the Nose−Hoover method 41 with Langevin piston pressure control 42 was employed to control fluctuations in the barostat. Harmonic restraints to the protein backbone (with a force constant of 1.0 kcal mol −1 Å −2 ) and side chains (with a force constant of 0.1 kcal mol −1 Å −2 ) were applied by means of the definition of a set of collective coordinates, 43 while the remaining system was kept free of restraints. Next, those harmonic restraints were fully released and a production step of 100 ns was run within the same (NPT) ensemble. An integration time step of 2 fs was used along the whole MD trajectory. All bonds and angles containing hydrogen atoms were constrained by the SHAKE algorithm. 44 Long-range electrostatic interactions were treated by the particle mesh Ewald (PME) method. 45 A cutoff of 12 Å was used to treat the nonbonding interactions.
For the Gly427Ala mutant the same protocol as described above was applied, except for the production step, for which a total of 50 ns were run.
2.5. QM/MM Calculations. The program Q-CHEM version 4.2 46,47 with the CHARMM interface 48 was used. The complete system derived from the MD simulations was trimmed to a 20 Å radius sphere centered on the AA molecule, containing over 4300 atoms. All residues and water molecules within a 15 Å radius sphere, including the complete AA molecule, were set as the active region (around 2030 atoms) of the optimization processes and allowed to move freely, whereas the remaining atoms were kept frozen. The QM region was described at the B3LYP-D3/6-31G(d,p) 49,50 level for all atoms. The MM region was described by the CHARMM36 39 force field. For the study of the hydrogen abstraction step (with a sextet multiplicity), the QM region consists of 55 atoms ( Figure 2a): 8 atoms of each imidazole ring of the three histidines, the side chain of Asn574, the terminal carboxylate group of Ile693, the pentadienyl group around C 10 of the lipid substrate, and the Fe(III)−OH − cofactor. A full electrostatic embedding scheme has been employed, and seven link atoms were used for the QM/MM boundary. No cutoffs were introduced for the nonbonding MM and QM/MM interactions.
For the oxygen insertion step (with a doublet multiplicity because the Fe atom is included in the MM region), a smaller QM region was employed, containing the oxygen molecule and the pentadienyl radical centered on C 10 , making a total of 18 QM atoms (Figure 2b). In this case, two link atoms were employed in the QM/MM boundary.
2.6. Free Energy Calculations. The free energy profile has been calculated using the umbrella sampling (US) method. For the hydrogen abstraction step, a total of 19 windows separated by 0.2 Å along the reaction coordinate were selected, comprising the entire range of the reaction coordinate. In the same way, 13 windows were selected to study the oxygen insertion step. For each window, an initial 0.2 ps phase of classic MD was performed, followed by 3 ps of QM(B3LYP-D3/6-31G(d,p))/MM dynamics. Each simulation was restrained at the corresponding particular window with a force constant of 250.0 kcal mol −1 Å −2 . An integration time step of 1 fs was used. The free energy profile was computed employing two different algorithms, the binless implementation of the weighted histogram analysis method (WHAM), 51 and a highprecision implementation of the dynamic histogram analysis method (DHAM) 52 to calculate both the free energy profile and the kinetic rate constant of the reaction. QM(B3LYP-D3/ 6-31G(d,p))/MM 2D free energy surfaces were also calculated (see below), involving 19, 29, and 48 umbrella windows (3 ps of production time for each one) for the hydrogen abstraction, the oxygen insertion to WT, and the oxygen insertion to the Gly427Ala mutant, respectively.

Hydrogen Abstraction
Step. 3.1.1. Molecular Dynamics Simulations. The starting structure for the MD simulations was derived from crystal structure 4QWT, subunit C, crystallized with the substrate in the active site. 18 In this structure the substrate is positioned in a tail-first orientation with the AA carboxylate end interacting with Arg182.
The first criterion to select a starting structure for the QM/ MM calculations in good agreement with the regioselectivity of the system is based on the distances of the three reactive bisallylic methylenes to the oxygen atom of the Fe(III)−OH − cofactor. These distances have been monitored along the MD simulation giving the results presented in Figures 3a and 3b. The average of the d(C 7 −OH), d(C 10 −OH), and d(C 13 −OH) distances (standard deviation in parentheses) is 7.57 (2.09) Å, 5.47 (1.33) Å, and 6.16 (1.50) Å, respectively; therefore, C 10 is on average closer to the catalytic iron than C 7 and C 13 , as would be expected for a 8-LOX reaction specificity. As for the evolution of the three distances along the MD simulation, during the first 15 ns, C 10 is clearly closer to the Fe(III)−OH − cofactor. After 15 ns, all three bisallylic methylenes are at greater distances, but d(C 10 −OH) remains as the shortest distance until 45 ns. Then an inversion of the distances order occurs and d(C 13 −OH) becomes the shortest distance until 85 ns. After that, d(C 10 −OH) becomes again the closest one. Thus, a fluctuation of the distances is produced during the MD production trajectory, but a good number of conformations where C 10 is in a better position with respect to the oxygen of the OH − are generated. Figures 4a and 4b show the first and last frames of the MD production step.
The RMSD of the protein heavy atoms with respect to the initial frame of the production step, and the RMSD of the AA heavy atoms with respect to the average structure are represented in Figure S2. A jump in protein RMSD is observed at 45 ns, which mainly corresponds to the fluctuation of the Nterminal domain. From this point, protein RMSD stabilizes. In the case of the AA RMSD, a larger value is observed at the beginning of the simulation, after which the RMSD also stabilizes.
3.1.2. QM/MM Calculations. Many of the MD configurations do not correspond to productive structures. Only those with a suitable (that is, enough short) pro-R H 10 −OH distance will be adequate to produce the hydrogen transfer step. Following this distances criterion, we have selected a particular suitable snapshot to start the QM/MM calculations and obtain the converged potential energy profile after repeated forward and backward processes, as explained below. The structure for the subsequent QM/MM calculations ( Figure S3) has been selected from the first 45 ns of the MD production step. Geometric parameters of interest of this structure are given in Table 1.
It is observed that the reactant structure fulfills the conditions of 8-lipoxygenation specificity in terms of hydrogen abstraction: the closest carbon atom to the oxygen of the Fe(III)−OH − cofactor is clearly C 10 ; pro-R H 10 is closer than C 10 , presenting thus a good orientation toward the reactive OH − to initiate the reactive process, and also, pro-S H 10 is further than pro-R H 10 (and C 10 as well), thus being not able to compete with pro-R H 10 for the reaction. This pro-R or pro-S nomenclature is assuming that the oxygen entrance is produced with antarafacial stereochemistry with respect to the face where hydrogen abstraction takes place, and taking into account that the main product comes from the oxygen addition at position C8; this pro-R H10 would lead to the oxidation product with S stereochemistry if the oxygen insertion took place at C 12 (see Figure S1 for details).
The optimized geometry of the reactant minimum has been used as the starting point to construct the potential energy profile along the reaction coordinate. The reaction coordinate, z, is defined as the difference between the distance of the breaking bond r 1 (C 10 −H 10 ) and the forming bond r 2 (H 10 −O), with H 10 being the pro-R hydrogen. To construct the potential energy profile, a series of optimizations of the mobile part of the system have been performed in the presence of harmonic restrictions on the reaction coordinate, which increases with a

ACS Catalysis
Research Article step size of 0.2 Å. The process has been repeated forward and backward until convergence to avoid hysteresis. 53 The resulting potential energy profiles are depicted in Figure 5.
Along the converged potential energy profile profile, reactants are located approximately at z = −1.3 Å, and products at z = 1.3 Å. For the transition state (TS) structure, z = 0.0 Å, which indicates that TS is highly symmetrical with respect to the movement of the transferred atom along the reaction coordinate. The potential energy barrier for pro-R H 10 abstraction is around 15.0 kcal/mol. This is a lower barrier than in previous studies of hydrogen abstraction processes in similar systems (19.6 kcal/mol in rabbit 15-LOX-1 24 and 18.0 kcal/ mol in human 15-LOX-2 26 for pro-S H 13 abstraction). This result is not surprising, since the catalytic constant for 8R-LOX is higher than for other LOX isoforms. The experimental k cat value for coral 8R-LOX is 210 ± 23 s −1 at room temperature, 18 whereas it is only 8.2 ± 0.25 and 1.5 ± 0.03 s −1 for rabbit 15-LOX-1 54 and human 15-LOX-2, 55 respectively. On the other hand, the reaction for 8R-LOX is exoergic by 10.0 kcal/mol.
3.1.3. Spin Densities. Atomic populations from spin density, which indicate the excess of α spin, have been evaluated during the reactive process for atoms C 8 , C 9 , C 10 , C 11 , and C 12 of the pentadienyl system, the abstracted hydrogen pro-R H 10 , and the catalytic iron (see Figure S4). Initially, the spin density on the mentioned AA atoms is zero, and as the reaction progresses, an increasing spin density appears on C 10 and delocalizes over C 10 , C 8 , and C 12 , going from 0 au to 0.5 au over each of these atoms. At the same time, the spin density on C 9 and C 11 changes from 0 to −0.2 au. The sum of spin densities over all these carbon atoms yields approximately 1.0 au, which corresponds to one unpaired electron delocalized over the pentadienyl radical system. The spin density on Fe decreases from 4.2 to 3.8 au, which corresponds to the transition of a Fe(III) sextet configuration to a Fe(II) quintet configuration. The spin density on pro-R H 10 remains virtually zero during the reactive process, indicating that this atom corresponds to a proton rather than a hydrogen atom bearing an electron. Furthermore, this atom is positively charged during the reaction process, with a charge changing from 0.15 to 0.38 au. Thus, this process corresponds to a proton coupled electron transfer (PCET) process, 56−60 in which the proton and the electron are transferred to different acceptors: the proton is transferred to the OH group oxygen to produce water, whereas the electron is transferred from the C 9 −C 13 pentadiene group of AA to the Fe(III) to form Fe(II).
3.1.4. Free Energy Calculations. To account for the thermal and entropic contributions, the free energy profiles were calculated using the umbrella sampling (US) protocol. Biased QM/MM molecular dynamics simulations were run using the same QM region as previously defined (see Figure 2). We used the recently developed DHAM method 52 to calculate both the free energies as well as the kinetic rates corresponding to the catalytic reaction. We used a new high precision Matlab

ACS Catalysis
Research Article implementation of the DHAM method 52 to allow us to numerically obtain accurate free energy profiles for the reaction with a high free energy barrier where the largest eigenvalues of the Markov matrix are almost degenerate. The histograms of the different US windows are represented in Figure S5, where it is observed that there is good overlap among the windows. The Markov matrices corresponding to the unbiased transition probabilities of the system have been calculated at ten different lag times (1, 2, 3, ..., 10 fs) and three different numbers of bins (600, 800, and 1000), and their spectral decomposition has been determined. The slowest relaxation rate corresponding to the second largest eigenvalue provided an average rate constant of 35 s −1 with a considerable standard deviation of 30 s −1 , calculated using the different conditions (lag times and number of bins) above. The free energy barrier using the DHAM method was found to be 16.0 ± 0.2 kcal/mol (Figure 6a), very similar to the potential energy barrier ( Figure  5). The DHAM calculated free energy compares well with the phenomenological free energy barrier of 14.4 kcal/mol at T = 300 K (k cat = 210 s −1 ). The unidimensional free energy profile using the traditional WHAM method has also been calculated for comparison with DHAM ( Figure 6a). A free energy barrier of 17.7 kcal/mol is obtained in this case. In addition, a twodimensional free energy surface as a function of the coordinates r 1 and r 2 , the two components of the reaction coordinate z, has been calculated using the WHAM method (see Figure 6b). All these free energy methods provide roughly identical results in this case.
At this point, it can be noted that, according to canonical variational transition-state theory, 61,62 the free energy barrier (ΔG CVT ) is related to the rate constant by means of the equation (without tunneling) where k is the rate constant, T is the temperature, k B is Boltzmann's constant, h is Planck's constant, R is the gas constant, and Γ is the quassiclassical transmission factor that corrects the rate constant for dynamical recrossing through the dividing surface corresponding to the canonical variational transition (CVT) state (that is, where the generalized transition-state free energy of activation is maximum). The frequency factor k B T/h equals 6.3 × 10 12 s −1 at T = 300 K, but the calculation of Γ is not easy. The advantage of using the DHAM method is that, at least at an approximate level, the barrier-crossing times can be estimated directly from the global analysis of local umbrella sampling trajectories, so that the effect of the dynamical recrossing is incorporated. Thus, the entire pre-exponential factor in eq 1 can be obtained substituting in that equation the rate constant and the free energy barrier calculated by the DHAM method, leading to an appropriate value of approximately 14.2 × 10 12 s −1 at 300 K. This value does not indicate a significant contribution of dynamical recrossing at the CVT dividing surface in this case. Such a conclusion is in agreement with the results obtained in a previous study carried out in our research group on the Habstraction process in 15-LOX-2. 26

Oxygen Insertion
Step. 3.2.1. Initial Position for Oxygen Attack. To start the study of the second reaction step, the oxygen addition to the AA pentadienyl moiety, the structure corresponding to the minimum of the product of the pro-R H 10 hydrogen abstraction step (see previous section), has been fully QM/MM optimized. The resulting structure consists of a radical delocalized over the C 8 −C 9 −C 10 −C 11 −C 12 pentadienyl system (see spin density in Figure S4). A procedure for positioning the oxygen molecule on this structure has been carried out. Several oxygen molecules were initially placed. These molecules have been positioned taking into account the Cartesian axes (Figure 7a) making a total of 18 different starting positions for the oxygen molecules (thus generating 18 initial positions for an oxygen molecule).

ACS Catalysis
Research Article Figures 7b and 7c show the initial poses of the oxygen molecules and the AA substrate. The oxygen molecules have been placed pointing to the target carbon atom with an initial distance from the closest oxygen atom to the carbon atom of 3.5 Å. Two different sets of structures were generated, those prepared for the oxygen attack at C 8 , and those corresponding to the attack at C 12 , making a total of 36 starting structures.
For each possible position generated for the O 2 molecule, a single point energy calculation has been performed in order to filter out the most stable ones. The results show that there are two positions that clearly exhibit higher stability (one corresponding to the placement around C 8 and the other to C 12 ), and it stands out that in both positions the oxygen molecule is similarly placed (see Figure S6), located above C 10 . In principle, oxygen addition to both C 8 and C 12 could be feasible. Thus, these two positions have been selected as the starting points of the reaction coordinate for oxygen attack at C 8 and C 12 , respectively. It is also interesting that the two most stable oxygen poses have a relative orientation with respect to the Fe(III)−OH − cofactor of antarafacial stereochemistry, which is in good agreement with the known experimental data. Thus, after the hydrogen abstraction step, the antarafacial oxygen addition would produce the 8-H(p)ETE oxidation product with R stereochemistry (8R-H(p)ETE), or the 12-H(p)ETE oxidation product with S stereochemistry (12S-H(p)ETE). These single point energy calculations were performed with the larger QM part (Figure 2a) also adding the O 2 molecule (with a sextet multiplicity).

Potential Energy Profiles for the Oxygen Insertion
Step. Following the same procedure as for the hydrogen abstraction step, the potential energy profiles for the oxygen attack have been constructed. The reaction coordinate for the second reaction step, z 2 , has been defined as the distance between the closest oxygen atom of the incoming O 2 molecule and the target carbon atom (C 8 or C 12 ) of the pentadienyl radical system. This coordinate has been scanned forward and backward until convergence to avoid hysteresis. 53 A step size of 0.2 Å has been defined. The potential energy profiles for the insertion at C 8 and C 12 are shown in Figures 8a and 8b, respectively. As the oxygen attack process only involves the oxygen molecule and the AA radical substrate, a small QM region containing only the pentadienyl moiety and the oxygen molecule (18 atoms), and excluding the Fe coordination environment, has been selected (Figure 2b). The QM region has been described at the B3LYP-D3/6-31G(d,p) level of theory.
For oxygen attack at C 8 , reactants lay approximately at z 2 = 3.0 Å, products at z 2 = 1.5 Å, and TS at z 2 = 2.2 Å. The potential energy barrier for this process is approximately 2.5 kcal/mol, and the reaction is exoergic by 3.6 kcal/mol. In the case of the oxygen addition at C 12 , reactants are located at z 2 = 2.9 Å, products at z 2 = 1.5 Å, and TS at z 2 = 2.1 Å. In this case, the potential energy barrier is 3.6 kcal/mol, and the reaction is exoergic by 1.6 kcal/mol. Thus, there is a difference of 1.1 kcal/ mol in potential energy barriers, with the insertion at C 8 being the most favorable process, which is in agreement with the experimental data. Atomic populations from spin density for the oxygen insertion step are also provided (see Figure S7).
In both cases, spin density on C 8 and C 12 increased from −0.5 to 0 au and that on the oxygen atoms decreased from 1 au in both atoms to 0.7 au in the case of O B and to 0.3 au in O A , which is the oxygen participating in the peroxide bond. Therefore, at the end of the process, the unpaired electron initially delocalized over the pentadienyl system delocalizes on the two oxygen atoms of the −OO • , as the sum of their spin densities yields 1 au (0.3 + 0.7).
3.2.3. Free Energy Profiles. Following the same procedure as for the hydrogen abstraction step, the free energy profiles have been computed for the oxygen insertion step using the US method. A total of 39 ps of production time distributed in 13 umbrella windows were run by QM/MM MD simulations at the B3LYP-D3/6-31G(d,p) level and employing the small QM region (Figure 2b) defined for the reaction coordinate of the oxygen insertion step. To compute the free energy profiles, both DHAM and WHAM methods were used. The histograms of the US windows are represented in Figure S8, showing a good overlap.
All the 1D free energy curves have been fitted to a ninth degree polynomial function, and the reported free energy barriers have been calculated from the maximum (TS region) and minimum of the reactant valley of the fitting curve (Figures  9a and 9b). The free energy barriers (Table 2) resulted in 5.1 and 5.2 kcal/mol for oxygen insertion at C 8 and C 12 , respectively, by the WHAM method, but they turned out to be 5.0 and 5.7 kcal/mol for attack at C 8 and C 12 , respectively, by the DHAM method. The difference of 0.7 kcal/mol (DHAM) is slightly smaller than the corresponding difference between the potential energy barriers, but is again in agreement Figure 8. Potential energy profiles for the oxygen insertion at C 8 (a) and C 12 (b) in WT as a function of the corresponding reaction coordinate z 2 (the distance from the closest oxygen atom up to C 8 or C 12 , respectively).

ACS Catalysis
Research Article with the experimental regio-and stereospecificities of the second reaction step.
The two 1D free energy curves have been calculated and fitted independently from each other. Since the precise determination of the corresponding free energy minima is difficult, both free energy barriers involve some uncertainties. Then, we can wonder if the small difference between them is significant enough to indicate which of the carbon atoms will be attacked by molecular oxygen. With the aim that the oxygen insertions have the same initial starting point, we performed a joint analysis of the free energy calculations to ensure that the two possible reaction outcomes use the same reference zero free energy. For this, an additional 2D free energy contour plot was obtained analyzing the data together for the two reactions (oxygen attack to C 8 and C 12 ) using WHAM, shown in the same graph (Figure 9c). It is worth noting here that transition state theory is not directly related to a 2D free energy surface, but to an 1D free energy curve. However, although rate constants cannot be calculated from a 2D free energy surface, these 2D surfaces permit us to understand the free energy landscape and shed light on the reaction mechanism. In this case, the 2D results show that the preferred reaction outcome at C 8 has a more direct pathway (entropically favored), as well as a lower free energy barrier (5 kcal/mol versus 6 kcal/mol for oxygen insertion at C 12 ) and it is also more exoergic.
According to this contour plot, the reaction through C 8 is again favored both thermodynamically and kinetically.
3.2.4. Discussion. During the approach of the oxygen molecule toward the target carbon atoms, protein residues close to the O 2 molecule and the carbon atoms have been analyzed. The Leu385 residue, which is close to C 8 but even closer to C 12 , seems to hinder the path of the oxygen toward C 12 , favoring the attack at C 8 . Leu385 is spatially close to Gly427, the residue proposed to control the stereochemistry of the reaction. Both residues are placed over the plane of the pentadienyl system. There are no protein atoms closer than 3 Å from C 8 , but Leu385 is closer than this distance from C 12 . A van der Waals surface representation of the protein residues close to the pentadienyl system shows that there is a channel in the vicinity of C 8 (Figures 10a, 10b, and 10c). This channel can also be observed by the tunnel computed by Caver. 35 This channel, which is enclosed by Gly427 and Leu385, connects C 8 with the region where the most stable initial oxygen positions were obtained.
As can be seen in Figure 10, the channel proposed here leads to C 8 rather than to C 12 , and this could explain the regioselectivity of the oxygen insertion step. Also, the position of Leu385 seems to be hindering the reaction path toward C 12 .
3.3. Study of the Gly427Ala Mutant. 3.3.1. MD Simulations. The mutant species Gly427Ala was prepared from the optimized structure of the hydrogen abstraction product of the WT system. For the mutant structure, additional MD simulation runs (50 ns) have been performed in order to reaccommodate the AA radical in the new environment. To run the MD simulations, the complete system (protein + waterbox) was recovered from the trimmed system employed in the previous QM/MM calculations. The parameters for the radical AA were derived from a previous work 30 using the charges obtained for the radical product in the corresponding QM/MM calculation.
Protein and AA RMSD of heavy atoms, calculated with respect to the initial and average structures of the MD production trajectory, respectively, are represented in Figure   Figure 9. DHAM (blue) and WHAM (red) free energy profiles for the oxygen insertion step at C 8 (a) and C 12 (b) in WT as a function of the corresponding reaction coordinate z 2 , respectively. (c) 2D free energy surface as a function of the two reaction coordinates z 2 (Q C8 and Q C12 for the oxygen addition to C 8 and C 12 , respectively) obtained by the WHAM method. The cross indicates the location of the starting O 2 free energy minimum.

ACS Catalysis
Research Article S9. Protein RMSD increases up to 3.5 Å. AA RMSD stabilizes at 1.5 Å during all the production steps, indicating that the fluctuations of the radical inside the binding cavity are small. A structure from the MD simulation of the Gly427Ala 8R-LOX:AA(radical) complex has been selected to perform the second reaction step. During the MD simulation process, a redistribution of the residues in contact with the mutated Ala427 takes place, leading to a new configuration of the active site cavity. To find the initial position of the oxygen molecule, the same procedure as described above has been carried out. Thus, a total of 36 oxygen poses (18 at C 8 and 18 at C 12 ) have been generated. Single point energy calculations lead to the selection of the positions shown in Figure S10.
Like in the WT system, the most stable position to attack C 12 in the mutant presents the oxygen molecule in an antarafacial configuration with respect the Fe(II)−OH 2 cofactor. The second position, prepared to react with C 8 , presents an oxygen molecule almost coplanar to the pentadienyl radical moiety; however, as the oxygen approaches to C 8 in the reaction coordinate, it yields an antarafacial configuration. Thus, these two positions have been selected as the starting points of the reaction coordinate for oxygen attack at C 12 and C 8 , respectively.
3.3.2. Potential Energy Profiles. The reaction coordinate settings were the same as for the WT system. The potential energy profiles corresponding to the oxygen addition at C 8 and C 12 of the two most stable positions are represented in Figures  11a and 11b, respectively.
Oxygen insertion at C 8 occurs with a potential energy barrier of approximately 4.1 kcal/mol and is exoergic by 8.4 kcal/mol. Reactants are located at z 2 = 2.7 Å, products at z 2 = 1.5 Å, and TS at z 2 = 2.1 Å. In the case of C 12 , the potential energy barrier is 1.3 kcal/mol, and the reaction energy −7.1 kcal/mol. Reactants lie at z 2 = 2.7 Å, products at z 2 = 1.5 Å, and TS at z 2 = 2.1 Å. Therefore, a dramatic inversion of the reaction specificity has been accurately reproduced, with the oxygenation process at C 12 favored over C 8 by 2.8 kcal/mol, in good agreement with the experimental results.
3.3.3. Free Energy Profiles. The free energy profiles for the oxygen insertion step in the Gly427Ala 8R-LOX mutant were calculated with the same procedure as for the WT system. The US windows histograms are represented in Figure S11, and the free energy profiles are shown in Figure 12.
For the oxygen insertion at C 8 , the 1D free energy barrier results in 5.8 (WHAM) and 6.0 (DHAM) kcal/mol, whereas, for the oxygen addition at C 12 , the 1D free energy barrier turns out to be 5.1 (WHAM) and 4.7 (DHAM) kcal/mol (see Figures 12a and 12b and Table 2). Again, as in the case of the WT system, this difference in the free energy barrier for both processes (attack at C 8 and C 12 ) decreases with respect to the difference in the potential energy barriers (which is of 2.8 kcal/ mol, favoring C 12 ), being now of 1.3 kcal/mol (DHAM), favoring the attack at C 12 .
As for the WT case, an additional 2D free energy contour plot was obtained for the two reactions (oxygen attack to C 8 and C 12 ) in the same graph (Figure 12c) by the WHAM method. The 2D free energy barriers for the oxygen insertion to C 8 and C 12 are 11.6 and 9.5 kcal/mol, respectively. This contour plot confirms that the O 2 addition to C 12 is now favored both thermodynamically and kinetically. Then, the inversion of the reaction specificity in the mutant is also obtained with the free energy calculations. It should be noted at this point that the 1D free energy barriers do not necessarily match the ones extracted from the 2D free energy surfaces, because they imply averaging over one more degree of freedom.
3.3.4. Discussion. To analyze the details of the changes produced by the introduction of a methyl moiety by the mutation of Gly427 into an alanine, a representation of the protein surface around the AA substrate has been depicted in Figures 13a and 13b. This representation corresponds to the starting structure of the mutant prior to the MD simulation. In these figures the mutant presents the same geometry as the WT system studied in the previous section (except for the mutated residue), as it corresponds to the starting structure used to generate the Gly427Ala 8R-LOX mutant. It can be seen that the additional methyl group in Ala427 with respect to Gly427 now Figure 10. van der Waals surface representation of protein residues close to the AA pentadienyl system (at < 5 Å) in WT. In (a), the representation has been superimposed with the initial position of the oxygen molecules. In (b), residues Leu385 and Gly427 are depicted. C 8 and C 12 are highlighted by a green and violet sphere, respectively. (c) Tunnel representation of the oxygen access channel toward C 8 .

ACS Catalysis
Research Article blocks the channel that led to C 8 in the WT. With the present position of the residues, the O 2 access to either C 8 or C 12 would be hindered.
However, during the MD simulation, the active site cavity reorganizes as a consequence of the mutation. This reorganization leads to a relocation of the oxygen access channel of the WT system, so favoring the oxygen entrance toward C 12 (Figures 14a and 14b). Then, a new different channel appears which is enclosed by Leu385, Ala427, Ile437, and Val384. This channel is directed closer toward C 12 , which would explain that change in the reaction specificity. Also, the placement of Ala427 pushes Leu385 and affects its surroundings, as explained in the Introduction. In an experimental work by Newcomer and co-workers 17 it has been proposed that Leu431 rather than Leu385 is the residue involved in the configuration of the oxygen access channel. In their crystallographic structure, Leu431 is indeed in closer contact with Gly427 than its counterpart Leu385, and it is also closer to the AA substrate. However, our MD simulations show that Leu431 moves to a further position during the simulations, whereas Leu385 remains closer to the AA substrate ( Figure S12). Both residues are conserved (Table S1), so it seems that the presence of a Leucine in this 385 position could be key for the regio-and stereocontrol of oxygenation.

CONCLUSIONS
In this paper we have combined MD simulations, QM/MM calculations, and umbrella sampling free energy simulations to study the reaction mechanism of the hydrogen abstraction from arachidonic acid by WT coral 8R-LOX, and oxygen insertion in both WT and the in silico mutant Gly427Ala of 8RLOX. The structure of 8R-LOX used here was taken from a crystal structure containing the AA substrate inside the active site cavity, which exhibits a tail-first orientation with its carboxylate end interacting with Arg182.
The MD simulations in WT 8R-LOX show that C 10 is on average the closest carbon atom to the catalytic iron, which is in agreement with the regioselectivity of 8R-LOX. During our MD simulations, the substrate retains its tail-first orientation and its interaction with the Arg182 residue.
The potential energy profile for the hydrogen abstraction step of pro-R H 10 occurs with a barrier of 15 kcal/mol. This barrier height is lower than in other LOX isoforms, which is in excellent agreement with the experimental result of a higher catalytic rate constant. The free energy barrier calculated by the DHAM method is of 16.0 ± 0.2 kcal/mol, very close to the potential energy barrier, which suggests that the entropic and thermal contributions for the first step in this system are small.
The second reaction step, the oxygen insertion, has been studied for the two possible positions starting from the H 10 abstraction step: C 8 (which should be the most favored position in the WT 8R-LOX), and C 12 . The potential energy barriers for the two processes turn out to be 2.5 and 3.6 kcal/mol, respectively. Thus, the oxygen addition to C 8 is favored by 1.1 kcal/mol, which agrees with the experimental regioselectivity. The oxygen addition to C 8 leads to the product with R stereochemistry, whereas addition to C 12 would yield the S product. The corresponding free energy barriers confirm that the addition at C 8 is the preferred one in WT. Analyzing the configuration of the active site cavity reveals a channel that is located over the pentadienyl moiety of the AA substrate, and that leads to C 8 rather than to C 12 attack. The channel is enclosed by Gly427 and Leu385, and Leu385 hinders the oxygen approach to C 12 .
When Gly427 is mutated to alanine, the configuration of the channel changes, being now redirected toward C 12 : Ala427 pushes Leu385, altering the configuration of the channel and blocking the region over C 8 . The potential energy barriers for the Gly427Ala mutant turn out to be 4.1 kcal/mol for oxygen addition to C 8 , and 1.3 kcal/mol for the oxygen addition to C 12 , and the corresponding free energy barriers corroborate that the regio-and stereospecificities are inverted with respect to the WT. All those results provide a detailed molecular view of the grounds of the Ala-versus-Gly concept. 6 It is worth noting that Newcomer et al. have proposed that Leu431 was the residue involved in the oxygen access channel. 17 However, our results show that Leu431 moves away from the AA substrate, whereas Leu385 remains closer. Interestingly, both Leu385 and Leu431 are highly conserved leucines. On the other hand, in a very recent paper, 14 which appeared while this one was being written, Kuḧn and Scheerer have used the crystal structures of the WT and its Ala420Gly mutant to propose Leu378 (that corresponds to Leu385 in coral 8R-LOX) as one of the relevant residues that governs the specificity in Pseudomonas aeruginosa LOX, another enzyme that holds the Ala-versus-Gly concept. The crystallographic data indicate that the mutation opens a new putative oxygen access tunnel which was closed in the WT. The presence of such type Figure 11. Potential energy profiles for the oxygen insertion at C 8 (a) and C 12 (b) in the mutant Gly427Ala 8R-LOX as a function of the corresponding reaction coordinate z 2 (the distance from the closest oxygen atom up to C 8 or C 12 , respectively).

ACS Catalysis
Research Article of residues (bulky and hydrophobic), conserved among many LOX isoforms, could be key in the chemical control 63 for the regio-and stereoselectivities of the reaction mechanism. The strategic position of the different leucines (or isoleucines) in the active site sterically hinders the substrate regions in a suitable way, resulting in the production of the final hydroperoxidation product with the desired regio-and stereospecificities. Atomistic level insights including dynamical effects are crucial for the complete understanding of how lipoxygenases achieve their extreme specificity.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acscatal.7b00842. Figure 12. DHAM (blue) and WHAM (red) free energy profiles for the oxygen insertion step at C 8 (a) and C 12 (b) in the mutant Gly427Ala 8R-LOX. (c) 2D-free energy surface as a function of the two reaction coordinates z 2 (Q C8 and Q C12 for the oxygen addition to C 8 and C 12 , respectively) obtained by the WHAM method. The cross indicates the location of the starting O 2 free energy minimum. Figure 13. van der Waals surface representation of protein residues close to the AA pentadienyl system (at <5 Å) in the starting structure of mutant Gly427Ala 8R-LOX. In (a) the representation has been superimposed with the initial position of the oxygen molecules used in the WT. In (b) residues Leu385 and Ala427 are depicted. C 8 and C 12 are highlighted by a green and violet sphere, respectively. Figure 14. (a) Surface representation of protein residues close to the AA pentadienyl system (at <5 Å) in the selected MD snapshot of mutant Gly427Ala 8R-LOX used to study the oxygen insertion step. Residues Leu385 and the mutated Ala427 are depicted. C 8 and C 12 are highlighted by a green and violet sphere, respectively. (b) Tunnel representation of the oxygen access channel toward C 12 .
Scheme of the reaction mechanism, arachidonic acid structure, MD simulations in WT Coral 8R-LOX, hydrogen abstraction step in WT Coral 8R-LOX, oxygen insertion step in WT Coral 8R-LOX, conserved residues of the active site cavity of 8R-LOX, free energy calculations of the oxygen insertion step in WT Coral 8R-LOX, MD simulations in the mutant Gly427Ala 8R-LOX, oxygen insertion step in the mutant Gly427Ala 8R-LOX, free energy calculations of the oxygen insertion step in the mutant Gly427Ala 8R-LOX, Leu385, and Leu431 (PDF)