Theoretical Characterization of the Step-by-Step Mechanism of Conversion of Leukotriene A4 to Leukotriene B4 Catalysed by the Enzyme Leukotriene A4 Hydrolase

LTA4H is a bifunctional zinc metalloenzyme that converts leukotriene A4 (LTA4) into leukotriene B4 (LTB4), one of the most potent chemotactic agents involved in acute and chronic inflammatory diseases. In this reaction, LTA4H acts as an epoxide hydrolase with a unique and fascinating mechanism, which includes the stereoselective attachment of one water molecule to the carbon backbone of LTA4 several methylene units away from the epoxide moiety. By combining Molecular Dynamics simulations and Quantum Mechanics/Molecular Mechanics calculations, we obtained a very detailed molecular picture of the different consecutive steps of that mechanism. By means of a rather unusual 1,7-nucleophilic substitution through a clear SN1 mechanism, the epoxide opens and the triene moiety of the substrate twists in such a way that the bond C6-C7 adopts its cis (Z) configuration, thus exposing the R face of C12 to the addition of a water molecule hydrogen-bonded to ASP375. Thus, the two stereochemical features that are required for the bioactivity of LTB4 appear to be closely related. The noncovalent π-π stacking interactions between the triene moiety and two tyrosines (TYR267 and, especially, TYR378) that wrap the triene system along the whole reaction explain the preference for the cis configuration inside LTA4H.


Introduction
Nowadays, it is recognised that inflammatory-based human diseases represent the leading causes of present-day morbidity and mortality worldwide [1]. Infection of the host, tissue injury or surgical trauma trigger the release of proinflammatory lipid and protein mediators by cells in the site of challenge. These chemical mediators act as chemoattractants to recruit neutrophiles. This initial inflammatory response occurs to protect the host and produces the very well know signs of acute inflammation (redness, heat, swelling and pain) [2]. After this onset phase, the inflammatory reaction must be resolved to prevent the inflammation from spreading. The resolution phase is regulated by a number of specialised proresolving lipid mediators (produced by human cells called macrophages). If resolution does not work well, the acute inflammation evolves to chronic inflammation [3][4][5][6][7]. Chronic inflammatory diseases are the most significant cause of death in the world today [1].
LTA4H is a monomeric bifunctional zinc metalloenzyme. It is detected in almost all mammalian cells, and it has both aminopeptidase activity and epoxide hydrolase activity. In this paper, we will focus on the second one. This epoxide hydrolase reaction is interesting because involves the stereoselective introduction of one water molecule to the carbon backbone (C12) of LTA4 in a position several methylene units away from the epoxy moiety (carbons C5 to C6) [15]. A high-resolution crystal structure of human LTA4H in the complex with the competitive inhibitor bestatin revealed that the protein includes three domains, an N-terminal domain, a Zn-containing catalytic domain and an α-helical C-terminal domain. The active site is placed in a deep cleft in between the domains, where the Zn 2+ is coordinated to HIS295, HIS299 and one carboxylic oxygen of GLU318 [16,17]. Because of the LTA4 instability (half-life of about 10 s at neutral pH) high-resolution crystal structures of LTA4H complexed to LTA4 were not determined until recently by Haeggström and coworkers [18]. They described six different structures of human LTA4H from five distinct crystal forms, finding two conformational states of the enzyme and several conformations of the substrate LTA4, and showing that LTA4H undergoes domain movements. From this structural information and previous site-directed mutagenesis experiments [15,19,20], Haeggström and coworkers [18] propose a mechanism for the epoxide hydrolase reaction. In short, the C1 carboxylate group of LTA4 anchors to ARG563, while Zn 2+ and TYR383 are coordinated to the epoxide oxygen. A catalytic water molecule (WAT1), polarised by Zn 2+ and GLU271, transfers a proton to produce an SN1 acid-induced opening of the epoxide ring. The positive charge of the resulting carbonium ion is delocalised over the triene system, with the C6-C7 bond in a pro-cis configuration stabilised by the L-shape of the tight pocket. Finally, a second water molecule, activated by ASP375, adds to C12 to generate an LTA 4 H is a monomeric bifunctional zinc metalloenzyme. It is detected in almost all mammalian cells, and it has both aminopeptidase activity and epoxide hydrolase activity. In this paper, we will focus on the second one. This epoxide hydrolase reaction is interesting because involves the stereoselective introduction of one water molecule to the carbon backbone (C 12 ) of LTA 4 in a position several methylene units away from the epoxy moiety (carbons C 5 to C 6 ) [15].
A high-resolution crystal structure of human LTA 4 H in the complex with the competitive inhibitor bestatin revealed that the protein includes three domains, an N-terminal domain, a Zn-containing catalytic domain and an α-helical C-terminal domain. The active site is placed in a deep cleft in between the domains, where the Zn 2+ is coordinated to HIS295, HIS299 and one carboxylic oxygen of GLU318 [16,17]. Because of the LTA 4 instability (half-life of about 10 s at neutral pH) high-resolution crystal structures of LTA 4 H complexed to LTA 4 were not determined until recently by Haeggström and coworkers [18]. They described six different structures of human LTA 4 H from five distinct crystal forms, finding two conformational states of the enzyme and several conformations of the substrate LTA 4 , and showing that LTA 4 H undergoes domain movements. From this structural information and previous site-directed mutagenesis experiments [15,19,20], Haeggström and coworkers [18] propose a mechanism for the epoxide hydrolase reaction. In short, the C 1 carboxylate group of LTA 4 anchors to ARG563, while Zn 2+ and TYR383 are coordinated to the epoxide oxygen. A catalytic water molecule (WAT1), polarised by Zn 2+ and GLU271, transfers a proton to produce an S N 1 acid-induced opening of the epoxide ring. The positive charge of the resulting carbonium ion is delocalised over the triene system, with the C 6 -C 7 bond in a pro-cis configuration stabilised by the L-shape of the tight pocket. Finally, a second water molecule, activated by ASP375, adds to C 12 to generate an R-hydroxyl group. It has to be underlined that the 12R-hydroxyl group and the 6Z,8E,10E configuration of the conjugated triene system are required for the bioactivity of LTB 4 [15]. It seems that ASP375 is needed for the generation of the former [15], whereas TYR378 could be involved in the formation of the right configuration of the triene moiety [21,22].
To our knowledge, only a theoretical study of the epoxide hydrolase reaction exists so far [23]. However, this work was just restricted to the step corresponding to the epoxy ring opening, using the self-consistent-charge density-functional tight-binding (SCC-DFTB) theory [24][25][26] for the Quantum Mechanics calculations and the CHARMM force field for the molecular mechanics part [27]. Thus, considering only the first reaction step and using a low level of calculation, it provided a rather incomplete view of the whole reaction mechanism to obtain LTB 4 .
In this paper, we combine Molecular Dynamics (MD) simulations and Quantum Mechanics/Molecular Mechanics (QM/MM) calculations, using the B3LYP hybrid functional to describe the QM part, to disclose the molecular details of the complete step-by-step mechanism of conversion of leukotriene A 4 to leukotriene B 4 catalysed by the enzyme leukotriene A 4 hydrolase. In particular, we analyse the key factors that drive the 12R and the 6Z,8E,10E stereochemistries of the hydroxyl group and the triene moiety, respectively, that are essential for the bioactivity of LTB 4 .

QM/MM Structural Analysis
According to the crystallographic structures [18], the catalytic water molecule WAT1 is hydrogen-bonded to two Glutamate residues, GLU271 and GLU296. However, the positions of the corresponding hydrogen atoms are not a priori clear. That is, which of these Glutamate residues, GLU271 or GLU296, is the protonated one is not known. Thus, as explained below, two versions of the protonation state of the LTA 4 :LTA 4 H Michaelis complex have been prepared starting from (GLU271, protonated GLU296) or (protonated GLU271, GLU296), leaving all the other protonation states untouched. These two initial structures have been QM/MM-optimised, and the same minimum energy structure was obtained in both cases, showing that the protonated residue is GLU296, as pictured in Figure 2, where the QM/MM partition used is also shown. Thus, WAT1 coordinates to GLU271 through a donor water hydrogen bond, while it connects to GLU296 through an acceptor hydrogen bond. The remaining water hydrogen atom is the one interacting with the oxygen atom of the epoxide. A complete view of this QM/MM-optimised LTA 4 :LTA 4 H Michaelis complex is pictured in Figure 3. A scheme of the main noncovalent interactions between the substrate LTA 4 and the enzyme LTA 4 H in this optimised Michaelis complex is shown in Figure 4. LTA 4 fits into the L-shaped hydrophobic cavity of the LTA 4 H. The substrate docks to the protein through the interaction between the polar head of the LTA 4 and ARG563 and LYS565 from the protein, which are located in the depth of the cavity. Thus, the substrate is placed in a head-first conformation. Moreover, this fitting conformation lets the epoxide oxygen move near the Zn environment, while it places C 12 near ASP375 (with WAT3 between them). The epoxide is not only stabilised by being in the second coordination sphere of the Zn atom, but also thanks to two hydrogen bonds from TYR383 and WAT1, the catalytic water which is also coordinated to Zn. With this conformation, TYR267 and TYR378 are the two residues closest to the triene moiety of the substrate. The main interactions between the substrate LTA 4 and the enzyme LTA 4 H in this optimised Michaelis complex are schematised in Figure 4. The role of all these residues is to limit the available space and give a shape to the cavity. TYR267 and TYR378 are key residues of the selectivity of the reaction too (see below).
The Zn ion is the central atom of the catalytic domain of the protein. In good agreement with the crystallographic structures [18], in the optimised structure it is pentacoordinated to two HISs (HIS295, HIS299), GLU318, WAT1 and WAT2 with an octahedral geometry (distance RMSD compared to the ideal octahedron: 0.490 Å) [28], so a vacancy for the coordination is available ( Figure 5). The free O from GLU318 is linked to WAT2 through a hydrogen bond. Moreover, there is a second coordination sphere, which initially contains the epoxide of LTA 4 , GLU271 (hydrogen-bonded to WAT1 and WAT2) and GLU296 (linked to WAT1 through a hydrogen bond).

Molecular Dynamics Simulation
In the above section we have analysed the structure of the optimised LTA 4 :LTA 4 H Michaelis complex and we have identified some key points for the reaction mechanism. However, we might wonder if these features correspond just to this particular minimum energy structure, or if they hold throughout the dynamic fluctuations of the enzyme system. To unravel this point, we analysed 100 ns extracted from a Molecular Dynamics simulation of that Michaelis complex carried out following the protocol explained below. The QM/MM-optimised structure described in Section 2.1. was taken as starting point for the MD simulation.
Firstly, we selected the snapshots along the Molecular Dynamics simulation where there is a water molecule with a distance ASP375-WAT up to 3 Å and a distance C 12 -WAT smaller than 5 Å at once. In Figure 6, we illustrate where the water molecule appears in those selected snapshots as a function of the distances ASP375-WAT and C 12 -WAT. This water molecule is the one that appears between ASP375 and C 12 . Along the simulation there is exchange between the water molecules, in such a way that the water molecule that is represented in this figure is not always the same. In the optimised QM/MM structure, this is WAT3. Almost all the snapshots (96% of the generated snapshots) include a bridging water molecule closer than 3 Å and 5 Å to ASP375 and C 12, respectively. Even as many as 58% of the generated snapshots involve a bridging water molecule closer than 3 Å and 4.5 Å to ASP375 and C 12, respectively. It is clear that in most configurations a water molecule, activated by ASP375, appears to be ready to attack C 12 , thus forming a hydroxyl group.   On the other hand, the evolution of distances shown in Figure 7 confirms that the epoxide of LTA 4 maintains two hydrogen bonds with TYR383 and WAT1, and that it remains in the second coordination sphere of Zn ion along the Molecular Dynamics simulation.
Finally, we can see in Figure 8 that the C 6 -C 7 bond of LTA 4 (the single bond next to the epoxide) keeps a rather pro-cis configuration along the Molecular Dynamics simulation. A total of 91% of the snapshots correspond to a pro-cis configuration |dihedral angle H6 − C6 − C7 − H7| ≤ 90 • , whereas only 9% are associated with a pro-trans configuration |dihedral angle H6 − C6 − C7 − H7| ≥ 90 • . The evolution of this dihedral angle along the complete enzyme reaction is a key point for the formation of LTB 4 , as will be explained below.

QM/MM Reaction Mechanism Calculations
The conversion of LTA 4 to LTB 4 is unique because the water molecule is proposed [15,17] to be introduced stereoespecifically at a site (C 12 ) six methylene units away from the epoxide moiety (C 5 to C 6 ), in what would be a 1,7-nucleophilic substitution. Although it is believed [17] to proceed by means of an S N 1 mechanism, an S N 2 mechanism should also be considered as possible. To unravel this point, a two-dimensional potential energy surface was built (see Figure 9) using a reaction coordinate r 1 c for the epoxide ring opening (r 1 c = d(C 6 (LTA 4 ) − O epox (LTA 4 )), that is, the epoxide breaking bond length, and a reaction coordinate r 2 c for the water addition to C 12 (r 2 , that is, the difference between the water breaking bond length and the forming O-C bond length. The rest of the geometrical parameters in the active region have been fully optimised. R stands for the optimised structure (Figures 3-5). Using the Dijsktra algorithm [29] we have also defined in Figure 6 the minimum energy reaction path (white points) on that surface and the minimum (R, INT1 and INT2) and maximum energy points along it (TS1 and TS2). LTA4 fits into the L-shaped hydrophobic cavity of the LTA4H. The substrate docks to the protein through the interaction between the polar head of the LTA4 and ARG563 and LYS565 from the protein, which are located in the depth of the cavity. Thus, the substrate is placed in a head-first conformation. Moreover, this fitting conformation lets the epoxide oxygen move near the Zn environment, while it places C12 near ASP375 (with WAT3 between them). The epoxide is not only stabilised by being in the second coordination sphere of the Zn atom, but also thanks to two hydrogen bonds from TYR383 and WAT1, the catalytic water which is also coordinated to Zn. With this conformation, TYR267 and TYR378 are the two residues closest to the triene moiety of the substrate. The main interactions between the substrate LTA4 and the enzyme LTA4H in this optimised Michaelis complex are schematised in Figure 4. The role of all these residues is to limit the available space and give a shape to the cavity. TYR267 and TYR378 are key residues of the selectivity of the reaction too (see below). As shown in Figure 9, the reaction takes place in two well-separated steps. Firstly, the C 6 − O epox bond progressively breaks, thus opening the epoxide ring. In the intermediate INT1, the epoxide is already open (the distance C 6 − Oepox is 2.3 Å), but the value of the reaction coordinate r 2 c has not significantly changed yet. From here on, in the second step the distance C6 − Oepox remains quite invariant, while a clear evolution of the reaction coordinate r 2 c occurs, indicating the water (WAT3) addition to C 12 and the breakage of an O-H bond in this water molecule and finally reaching the intermediate INT2. Thus, the 1,7-nucleophilic substitution clearly takes place through an S N 1 mechanism. As a matter of fact, it is an extreme case of S N 1 reaction, because the minimum energy reaction paths corresponding to each step appear to be roughly orthogonal (each one parallel to one coordinate axis) in the two-dimensional potential energy surface shown in Figure 9. Any deviation of the mechanism towards the region near the diagonal of that energy surface (that is, any approach to the S N 2 mechanism) implies the penetration in high-energy regions, which turns out to be forbidden. An NBO charge analysis in INT1 gives a charge of +0.71 a.u. delocalised over the triene system (from C 6 to C 12 ) of the substrate, which confirms the existence of the carbocation required in the intermediate of an S N 1 mechanism.
The Zn ion is the central atom of the catalytic domain of the protein. In good agreement with the crystallographic structures [18], in the optimised structure it is pentacoordinated to two HISs (HIS295, HIS299), GLU318, WAT1 and WAT2 with an octahedral geometry (distance RMSD compared to the ideal octahedron: 0.490 Å) [28], so a vacancy for the coordination is available ( Figure 5). The free O from GLU318 is linked to WAT2 through a hydrogen bond. Moreover, there is a second coordination sphere, which initially contains the epoxide of LTA4, GLU271 (hydrogen-bonded to WAT1 and WAT2) and GLU296 (linked to WAT1 through a hydrogen bond).

Molecular Dynamics Simulation
In the above section we have analysed the structure of the optimised LTA4:LTA4H Michaelis complex and we have identified some key points for the reaction mechanism. However, we might wonder if these features correspond just to this particular minimum energy structure, or if they hold throughout the dynamic fluctuations of the enzyme system. To unravel this point, we analysed 100 ns extracted from a Molecular Dynamics simulation of that Michaelis complex carried out following the protocol explained below. The QM/MM-optimised structure described in Section 2.1. was taken as starting point for the MD simulation.
Firstly, we selected the snapshots along the Molecular Dynamics simulation where there is a water molecule with a distance ASP375-WAT up to 3 Å and a distance C12-WAT smaller than 5 Å at once. In Figure 6, we illustrate where the water molecule appears in those selected snapshots as a function of the distances ASP375-WAT and C12-WAT. This water molecule is the one that appears between ASP375 and C12. Along the simulation there is exchange between the water molecules, in such a way that the water molecule that is represented in this figure is not always the same. In the optimised QM/MM structure, this is WAT3. Almost all the snapshots (96% of the generated snapshots) include a bridging water molecule closer than 3 Å and 5 Å to ASP375 and C12, respectively. Even as many as 58% of the generated snapshots involve a bridging water molecule closer than 3 Å and 4.5 Å to ASP375 and C12, respectively. It is clear that in most configurations a water molecule, activated by ASP375, appears to be ready to attack C12, thus forming a hydroxyl group.  On the other hand, the evolution of distances shown in Figure 7 confirms that the epoxide of LTA4 maintains two hydrogen bonds with TYR383 and WAT1, and that it remains in the second coordination sphere of Zn ion along the Molecular Dynamics simulation. placed between ASP375 and C12.
On the other hand, the evolution of distances shown in Figure 7 confirms that the epoxide of LTA4 maintains two hydrogen bonds with TYR383 and WAT1, and that it remains in the second coordination sphere of Zn ion along the Molecular Dynamics simulation. Finally, we can see in Figure 8 that the C6-C7 bond of LTA4 (the single bond next to the epoxide) keeps a rather pro-cis configuration along the Molecular Dynamics simulation. A total of 91% of the snapshots correspond to a pro-cis configuration | ℎ 6 6 7 7| 90°, whereas only 9% are associated with a protrans configuration | ℎ 6 6 7 7| 90°. The evolution of this dihedral angle along the complete enzyme reaction is a key point for the formation of LTB4, as will be explained below. whereas | ℎ | 90° indicates a pro-trans configuration (blue line).

QM/MM Reaction Mechanism Calculations
The conversion of LTA4 to LTB4 is unique because the water molecule is proposed [15,17] to be introduced stereoespecifically at a site (C12) six methylene units away from the epoxide moiety (C5 to C6), in what would be a 1,7-nucleophilic substitution. Although it is believed [17] to proceed by means of an SN1 mechanism, an SN2 mechanism should also be considered as possible. To unravel this point, a two-dimensional potential energy surface was built (see Figure 9) using a reaction coordinate for the epoxide ring opening ( = d(C6(LTA4)−Oepox(LTA4)), that is, the epoxide breaking bond length, and a reac-  As shown in Figure 9, the reaction takes place in two well-separated steps. Firstly, the C6 − Oepox bond progressively breaks, thus opening the epoxide ring. In the intermediate INT1, the epoxide is already open (the distance C6 − Oepox is 2.3 Å), but the value of the reaction coordinate has not significantly changed yet. From here on, in the second step the distance C6 − Oepox remains quite invariant, while a clear evolution of the reaction coordinate occurs, indicating the water (WAT3) addition to C12 and the breakage of an O-H bond in this water molecule and finally reaching the intermediate INT2. Thus, the 1,7-nucleophilic substitution clearly takes place through an SN1 mechanism. As a matter of fact, it is an extreme case of SN1 reaction, because the minimum energy reaction paths corresponding to each step appear to be roughly orthogonal (each one parallel to one coordinate axis) in the two-dimensional potential energy surface shown in Figure 9. Any deviation of the mechanism towards the region near the diagonal of that energy surface (that is, any approach to the SN2 mechanism) implies the penetration in high-energy regions, which turns out to be forbidden. An NBO charge analysis in INT1 gives a charge of +0.71 a.u. delocalised over the triene system (from C6 to C12) of the substrate, which confirms the existence of the carbocation required in the intermediate of an SN1 mechanism.
To our surprise, a structural analysis of the intermediate INT2 shows that the conjugated triene moiety has a 6E,8E,10E configuration that would not lead to LTB4 but to its stereoisomer 12S-6-trans-LTB4 (5S,12S-dihydroxy-6E,8E,10E,14Z-eicosatetraenoic acid). Following the evolution of the structures along the SN1 mechanism shown in Figure 9, we found that an unexpected rotation of the C6-C7 bond from the initial rather pro-cis configuration in LTA4 to a pro-trans configuration in INT1 occurs, which leads to the 6E inadequate configuration of C6-C7 in INT2. To analyse this fact, we built a two-dimensional potential energy surface (see Figure 10) as a function of the rotation of the C6-C7 bond ( = dihedral angle H6-C6-C7-H7) and the epoxide ring opening ( = d(C6(LTA4)−Oepox(LTA4)). Starting from the reactant R, the Dijsktra algorithm [29] localises two minimum energy reaction paths (grey and white points) on the surface shown in Figure 10, and the minimum (R, INT1trans and INT1cis) and maximum (TS1trans and TS1cis) energy points along them. One path involves the pro-trans (grey points) configuration (as a matter of fact, it corresponds to the first step of the SN1 mechanism shown in Figure 9). The other To our surprise, a structural analysis of the intermediate INT2 shows that the conjugated triene moiety has a 6E,8E,10E configuration that would not lead to LTB 4 but to its stereoisomer 12S-6-trans-LTB 4 (5S,12S-dihydroxy-6E,8E,10E,14Z-eicosatetraenoic acid). Following the evolution of the structures along the S N 1 mechanism shown in Figure 9, we found that an unexpected rotation of the C 6 -C 7 bond from the initial rather pro-cis configuration in LTA 4 to a pro-trans configuration in INT1 occurs, which leads to the 6E inadequate configuration of C 6 -C 7 in INT2. To analyse this fact, we built a two-dimensional potential energy surface (see Figure 10) as a function of the rotation of the C 6 -C 7 bond (r 3 c = dihedral angle H 6 -C 6 -C 7 -H 7 ) and the epoxide ring opening (r 1 c = d(C 6 (LTA 4 ) − O epox (LTA 4 )). Starting from the reactant R, the Dijsktra algorithm [29] localises two minimum energy reaction paths (grey and white points) on the surface shown in Figure 10, and the minimum (R, INT1trans and INT1cis) and maximum (TS1trans and TS1cis) energy points along them. One path involves the pro-trans (grey points) configuration (as a matter of fact, it corresponds to the first step of the S N 1 mechanism shown in Figure 9). The other one maintains a pro-cis configuration (white points) to reach an intermediate INT1cis with the adequate pro-cis structure. The two paths are competitive, with the pro-trans being the one that implies the lowest potential energy barrier. This explains why only the trans configuration was obtained in the two-dimensional potential energy surface shown in Figure 9, where the rotation of the C 6 -C 7 bond was not included to define the reaction coordinate. This result also agrees with our B3LYP/6-31G(d) calculation that shows that in the gas phase, 12S-6-trans-LTB 4 is 0.3 kcal/mol more stable than LTB 4 in terms of potential energy. However, the experimental result is that LTA 4 H converts LTA 4 into LTB 4 , that is, with a 6Z,8E,10E configuration for the conjugated triene moiety.
An especially interesting point is that the configuration of the C 6 -C 7 bond in INT1 determines the stereochemistry of the ingoing water into C 12 . In Figure 11, we have displayed the position of WAT3 with respect to the plane generated by the triene moiety (from C 6 to C 12 ) in the reactant LTA 4 (pro-cis, yellow disk) in INT1trans (pro-trans, brown disk) and in INT1cis (pro-cis, violet disk). It can be clearly seen that the face of attack of WAT3 in the pro-trans configuration is just the opposite than in the case of the pro-cis configuration in INT1, leading to the 12S-hydroxyl or the 12R-hydroxyl configurations, respectively. Thus, the two main stereochemical features that are required for the bioactivity of LTB 4 appear to be closely linked.
one maintains a pro-cis configuration (white points) to reach an intermediate INT1cis with the adequate pro-cis structure. The two paths are competitive, with the pro-trans being the one that implies the lowest potential energy barrier. This explains why only the trans configuration was obtained in the two-dimensional potential energy surface shown in Figure  9, where the rotation of the C6-C7 bond was not included to define the reaction coordinate. This result also agrees with our B3LYP/6-31G(d) calculation that shows that in the gas phase, 12S-6-trans-LTB4 is 0.3 kcal/mol more stable than LTB4 in terms of potential energy. However, the experimental result is that LTA4H converts LTA4 into LTB4, that is, with a 6Z,8E,10E configuration for the conjugated triene moiety. An especially interesting point is that the configuration of the C6-C7 bond in INT1 determines the stereochemistry of the ingoing water into C12. In Figure 11, we have displayed the position of WAT3 with respect to the plane generated by the triene moiety (from C6 to C12) in the reactant LTA4 (pro-cis, yellow disk) in INT1trans (pro-trans, brown disk) and in INT1cis (pro-cis, violet disk). It can be clearly seen that the face of attack of WAT3 in the pro-trans configuration is just the opposite than in the case of the pro-cis configuration in INT1, leading to the 12S-hydroxyl or the 12R-hydroxyl configurations, respectively. Thus, the two main stereochemical features that are required for the bioactivity of LTB4 appear to be closely linked. At this point, it is clear that a complete and comparative study of the pro-cis and the pro-trans reaction paths is needed to understand the mechanism of this enzyme reaction. To this aim, firstly we fully optimised the structures INT1cis, INT1trans, TS1cis and TS1trans, which appear in Figure 10  At this point, it is clear that a complete and comparative study of the pro-cis and the pro-trans reaction paths is needed to understand the mechanism of this enzyme reaction. To this aim, firstly we fully optimised the structures INT1cis, INT1trans, TS1cis and TS1trans, which appear in Figure 10, in order to locate the corresponding stationary structures (minima and transition-state structures) in the complete potential hypersurface of the reaction. Then, for each reaction path, we followed the reaction coordinates d

(O(WAT3)-H(WAT3))-d(O(WAT3)-C 12 (LTA 4 )) and d(O(WAT1)-H(WAT1))-d(H(WAT1)-
O epox ) to describe, respectively, the water addition (WAT3) to C 12 and the final protonation of the oxygen atom bonded to C 6 (the former oxygen atom of the epoxide) by WAT1 and to locate the corresponding stationary structures. The relative potential energies for all the stationary points we localised are shown in Figure 12. Figure 11. Position of WAT3 with respect to the plane generated by the triene moiety in the reactant LTA4 (pro-cis, yellow disk), in INT1trans (pro-trans, brown disk) and in INT1cis (pro-cis, violet disk).
At this point, it is clear that a complete and comparative study of the pro-cis and the pro-trans reaction paths is needed to understand the mechanism of this enzyme reaction. To this aim, firstly we fully optimised the structures INT1cis, INT1trans, TS1cis and TS1trans, which appear in Figure 10, in order to locate the corresponding stationary structures (minima and transition-state structures) in the complete potential hypersurface of the reaction. Then, for each reaction path, we followed the reaction coordinates d(O(WAT3)-H(WAT3))-d(O(WAT3)-C12(LTA4)) and d(O(WAT1)-H(WAT1))d(H(WAT1)-Oepox) to describe, respectively, the water addition (WAT3) to C12 and the final protonation of the oxygen atom bonded to C6 (the former oxygen atom of the epoxide) by WAT1 and to locate the corresponding stationary structures. The relative potential energies for all the stationary points we localised are shown in Figure 12. We focused firstly on the pro-cis reaction path. In Figure 13 we show a succession of pictures of the stationary structures along this path, intending to capture LTA 4 H in action. For the sake of clarity, in Figures S1-S3 we display representations of those stationary structures focused on the region where the corresponding step occurs in each case. We focused firstly on the pro-cis reaction path. In Figure 13 we show a succession of pictures of the stationary structures along this path, intending to capture LTA4H in action. For the sake of clarity, in Figures S1-S3 we display representations of those stationary structures focused on the region where the corresponding step occurs in each case. The evolution of the main geometrical parameters along the pro-cis reaction path leading to LTB4 is indicated in Table 1. As indicated above, the first step consists of the epoxide opening. The formation of a nascent negative charge in the oxygen of the epoxide is stabilised by the interaction of this oxygen with the Zn atom, the hydroxyl group of TYR383 and the hydrogen bond interaction with WAT2 and WAT1, although WAT1 The evolution of the main geometrical parameters along the pro-cis reaction path leading to LTB 4 is indicated in Table 1. As indicated above, the first step consists of the epoxide opening. The formation of a nascent negative charge in the oxygen of the epoxide is stabilised by the interaction of this oxygen with the Zn atom, the hydroxyl group of TYR383 and the hydrogen bond interaction with WAT2 and WAT1, although WAT1 moves away from INT1cis. Conversely, the Zn-O epox bond is significantly shortened from 3.39 Å in R to 1.99 Å at INT1cis, in such a way that epoxide completes the coordination sphere of Zn by occupying its free vacancy (see Figure 5). The epoxide is fully opened at INT1cis, while O epox -C 5 is stronger at INT1cis than in R. This step implies a potential energy barrier of 11.6 kcal/mol. As explained above, at INT1cis no significant changes have occurred with respect to WAT3 and ASP375 yet, thus indicating an S N 1 mechanism. Very interestingly, the dihedral angle H 6 -C 6 -C 7 -H 7 goes from −65.6 • in R (a value close to an intermediate situation of the bond C 6 -C 7 between pro-cis and pro-trans) to −15.9 • at TS1cis and −7.1 • at INT1cis, a clearly cis stereochemistry. As a result of the twist that the triene moiety must perform to adopt this cis stereochemistry, WAT3 becomes ready to attack the R face of C 12 (by the side opposite to the Zn atom, see Figure S2) at INT1cis. Then, WAT3 adds to C 12 (the distance O(WAT3)-C 12 evolves from 3.69 Å at INT1cis to 1.53 Å at INT2cis), while one of its hydrogen atoms is fully transferred to ASP375 through the corresponding hydrogen bond. This second step involves a potential energy barrier of 15.9 kcal/mol. Finally, in the third step, the final transfer of a proton from WAT1 to the oxygen atom bonded to C 6 (the former oxygen atom of the epoxide, O epox ) takes place, while the distance Zn-O epox clearly increases through a potential energy barrier of 14.4 kcal/mol, thus forming LTB 4 . Note that Haeggström and coworkers [18] have proposed that the catalytic water WAT1 transfers a proton to promote an S N 1 acid-induced opening of the epoxide ring. However, we have found here that this proton donation is not the first step of the reaction, but the last one. That is, when the epoxide ring opens, the nascent negative charge in its oxygen atom is quite stabilised through the interaction with atoms surrounding it, in such a way that a previous proton transfer from WAT1 is not needed. Table 1. Most relevant distances (in Å) and dihedral angle (in degrees) defining the E/Z character of the C 6 -C 7 bond for the stationary structures of the complete mechanism of conversion of LTA 4 to LTB 4 (pro-cis reaction path) catalysed by LTA 4 H. R  TS1  INT1  TS2  INT2  TS3  Let us describe the pro-trans reaction path now. In Figure 14, we show a succession of pictures of the stationary structures along this path. For the sake of clarity, in Figures S4-S6, we display representations of those stationary structures focused on the region where the corresponding step occurs in each case. Let us describe the pro-trans reaction path now. In Figure 14, we show a succession of pictures of the stationary structures along this path. For the sake of clarity, in Figures  S4-S6, we display representations of those stationary structures focused on the region where the corresponding step occurs in each case. The evolution of the main geometrical parameters along the pro-trans reaction path is indicated in Table 2. Analogously to the case of the pro-cis reaction path, in the first step of the pro-trans reaction path, the oxygen of the epoxide is stabilised by the interaction with the Zn atom, the hydroxyl group of TYR383 and the hydrogen bond interaction with WAT2 and WAT1, although WAT1 also moves away at INT1trans. Conversely, the Zn- The evolution of the main geometrical parameters along the pro-trans reaction path is indicated in Table 2. Analogously to the case of the pro-cis reaction path, in the first step of the pro-trans reaction path, the oxygen of the epoxide is stabilised by the interaction with the Zn atom, the hydroxyl group of TYR383 and the hydrogen bond interaction with WAT2 and WAT1, although WAT1 also moves away at INT1trans. Conversely, the Zn-O epox bond is significantly shortened from 3.39 Å in R to 1.98 Å at INT1trans. The epoxide is fully opened at INT1trans of this S N 1 mechanism. This step implies a potential energy barrier of only 7.9 kcal/mol. The key stereochemical features of the complete reaction are determined in this first step, where the fate of the process is decided. Thus, the dihedral angle H 6 -C 6 -C 7 -H 7 is 196.9 • at TS1trans and 183.7 • at INT1trans, a clearly trans stereochemistry. This twist of the triene moiety takes place in the opposite direction to the twist that happens in the pro-cis reaction path. This way, WAT3 now becomes ready to attack the S face of C 12 (by the same side as the Zn atom, see Figure S5) at INT1trans. Then, WAT3 adds to C 12 (the distance O(WAT3)-C 12 evolves from 3.66 Å at INT1trans to 1.51 Å at INT2trans), while one of its hydrogen atoms is fully transferred to ASP375 through the corresponding hydrogen bond. This second step involves a potential energy barrier of 18.7 kcal/mol. Finally, in the third step, the final transfer of a proton from WAT1 to the oxygen atom bonded to C 6 takes place, while the distance Zn-O epox clearly increases through a potential energy barrier of 17.8 kcal/mol, thus forming 12S-6-trans-LTB 4 , the stereoisomer of LTB 4 . Table 2. Most relevant distances (in Å) and dihedral angle (in degrees) defining the E/Z character of the C 6 -C 7 bond for the stationary structures of the pro-trans reaction path catalysed by LTA 4 H. R  TS1  INT1  TS2  INT2  TS3  A final point to discuss is why LTA 4 H converts LTA 4 into LTB 4 instead of the stereoisomer 12S-6-trans-LTB 4 , which is more stable in the gas phase. In other words, why the pro-cis reaction path is the dominant one instead of the pro-trans reaction path. As seen in Figure 12, as for the epoxide opening the trans path is more favourable than the cis path. However, from here on, the cis energy profile appears clearly below the trans energy profile. In particular, taking into account the complete mechanism, the higher energy transition state structure of the trans path (18.7 kcal/mol) and the trans final product (10.3 kcal/mol) are significantly above the corresponding values of the cis path (15.9 kcal/mol and 8.5 kcal/mol, respectively), what means that the inside of the LTA 4 H formation of LTB 4 is more favourable than the formation of 12S-6-trans-LTB 4 both kinetically and thermodynamically. This energy crossing can be explained by observing some of the most important noncovalent interactions that take place between the substrate and the enzyme. TYR267 and TYR378 are the two residues closest to the triene system, which is practically wrapped by these two tyrosines along the reaction paths (see the sequence of structures shown in Figures 13 and 14). The key point here is that these two tyrosines and the triene moiety can interact among them through noncovalent π-π stacking interactions that are different depending on the reaction path (see Figure 15 for the angles between the respective planes). The angle of the π-π interaction between TYR267 and TYR378 is kept quite invariant (between 50 • and 60 • ) throughout the complete reaction for both reaction paths, and it has not been depicted. It has to be recalled here that the closer the planes to the face-to-face orientation, the more stabilising the interaction. Thus, TYR267 slightly favours the pro-trans reaction path up to INT1 (the first step), but the pro-cis reaction path in the second and third steps. Much more relevant is the role of TY378. As explained above, during the epoxide opening the triene moiety twists in two opposite directions depending on which reaction path it takes, pro-cis or pro-trans. As a result of this twist, the π-π stacking interaction with TY378 becomes clearly face-to-face (pro-trans) or clearly edge-to-face (pro-cis). At INT1trans (see Figure 11) the planes of the triene moiety and TY378 are roughly parallel. This way, TYR378 favours the formation of a broken epoxide with a trans C 6 -C 7 bond. However, along the second and third steps, edge-to-face interaction is kept for the pro-cis reaction path, while the face-to-face interaction for the pro-trans reaction path becomes significantly broken, thus destabilising the formation of 12S-6-trans-LTB 4 and leading to LTB 4 as a product of the reaction catalysed by LTA 4 H. Indeed, those two tyrosines do not exist in the gas phase, where the 12S-6-trans-LTB 4 is more stable than LTB 4 . Our theoretical result agrees with experimental data [21,22] showing that mutations of TY378 lead not only to LTB 4 but also to products with a trans stereochemistry in the C 6 -C 7 bond.  Figure 15. Angles between the plane of the triene moiety of the substrate and the planes of the aromatic rings of TYR267 and TYR378 for the different stationary structures located along the pro-cis and the pro-trans reaction paths. The closest plane to the whole number of C and H atoms of the triene system was taken as the plane of this triene moiety.

Protein Setup
Numerous crystallographic structures [18] have been reported for human LTA4H, some of them complexing a wide variety of inhibitors and a few of them containing the LTA4 substrate but including one or more mutations. Here, a single-mutated version of human LTA4H (ASP375ASN) complexed with LTA4 was used (PDB code: 5NI6) [18]. The mutation was reverted to the WT enzyme. The protein was protonated at pH 7.0 using PropKa3.0 [30] through a web interface (www.playmolecule.org, access on 1 July 2019). Two versions of the protonation state were prepared: (GLU271, protonated GLU296) and (protonated GLU271, GLU296), leaving all the other protonation states untouched. A box of TIP3P [31] water molecules was built around the protein for solvation. The box was built up considering 10 Å from the outermost protein's atom in each direction of the space.

QM/MM Calculations
The solvation model was cropped, so only the water molecules inside a 17 Å radius sphere around whatever atom of the substrate were kept for the QM/MM calculations. The modular software package ChemShell [32,33] was used as the interface between the QM and the MM calculations. TURBOMOLE [34] was used for the density functional the- Figure 15. Angles between the plane of the triene moiety of the substrate and the planes of the aromatic rings of TYR267 and TYR378 for the different stationary structures located along the pro-cis and the pro-trans reaction paths. The closest plane to the whole number of C and H atoms of the triene system was taken as the plane of this triene moiety.

Protein Setup
Numerous crystallographic structures [18] have been reported for human LTA 4 H, some of them complexing a wide variety of inhibitors and a few of them containing the LTA 4 substrate but including one or more mutations. Here, a single-mutated version of human LTA 4 H (ASP375ASN) complexed with LTA 4 was used (PDB code: 5NI6) [18]. The mutation was reverted to the WT enzyme. The protein was protonated at pH 7.0 using PropKa3.0 [30] through a web interface (www.playmolecule.org, accessed on 1 July 2019). Two versions of the protonation state were prepared: (GLU271, protonated GLU296) and (protonated GLU271, GLU296), leaving all the other protonation states untouched. A box of TIP3P [31] water molecules was built around the protein for solvation. The box was built up considering 10 Å from the outermost protein's atom in each direction of the space.

QM/MM Calculations
The solvation model was cropped, so only the water molecules inside a 17 Å radius sphere around whatever atom of the substrate were kept for the QM/MM calculations. The modular software package ChemShell [32,33] was used as the interface between the QM and the MM calculations. TURBOMOLE [34] was used for the density functional theory (DFT) calculations, while AMBER [35] force fields were employed for the MM calculations using the DL_POLY [36] module in ChemShell. AMBER's ff14SB [37] force field was used for generating the protein MM parameters. GAFF2 [38] force field was used for parameterising the substrate (LTA 4 ), whose RESP charges [39] were calculated using Gaussian09 [40]. An electrostatic embedding scheme [41] was employed in order to treat the QM mm interactions. Additionally, a link atom scheme was used to treat the QM/MM boundary by using the charge shift model [42]. No cut-offs were introduced for the nonbonding QM/MM and MM interactions [43].
The active region of the system was built by selecting all atoms within a radius of 15 Å around the oxygen atom of the epoxide group of the substrate. These atoms were allowed to move freely (almost 3000 atoms), while the nonactive ones were frozen (almost 8000 atoms). For the QM part of the system (see Figure 2), 86 atoms were selected involving the Zn and its immediate environment (HIS295, HIS299, GLU318, WAT1, WAT2), GLU271 and GLU296, from C 4 to C 13 of the substrate (thus including the epoxide and the triene group and ASP375 and its nearest water molecule (WAT3). From this selection, eight link atoms were added: six between the corresponding bonds C(MM)-C(QM) atoms of the six residues and two bonded to the aliphatic carbon atoms of the lipid substrate (placed between C 3 -C 4 and C 13 -C 14 ). Moreover, a microiterative scheme [44] was applied for the energy optimisation to minima. All residues containing at least one atom in the QM region were added to the 'micro' region of the microiterative scheme, while all other atoms were kept in the 'macro' region.
The QM/MM energy optimisations to minima and scan calculations were performed using the HDLCopt (hybrid delocalised internal coordinate) coordinates scheme [45] and the limited-memory Broyden-Fletcher-Goldfard-Shanno (L-BFGS) algorithm [46,47]. For optimisations to minima, tolerance was set to 0.0045 Bohr, while for potential energy surface (PES) calculations it was set to 0.01 Bohr. PES calculations have at least one angle, bond or bond distance differences fixed in order to carry out the exploration. The transition-state searches were carried out employing the HDLC coordinates scheme and the partitioned rational function optimiser (P-RFO) [48,49] combined with the L-BFGS algorithm. P-RFO and the L-BFGS algorithm were used as implemented in the HDLCopt module and the DL_FIND geometry optimisation library [50] of ChemShell, respectively. A set of core atoms including only the 6 atoms directly implicated in the reaction was defined in order to ease the calculation of the Hessian along the optimisation. Frequencies of the optimised structures were calculated using the force module from ChemShell for the whole QM region. The nature of the stationary structures was confirmed by means of the analysis of the number of imaginary frequencies.
We used the Dijsktra algorithm [29] to obtain the minimum energy reaction path on the two-dimensional potential energy surfaces.
The QM region was described by the B3LYP hybrid functional [51]. The 6-31G(d) Pople basis set [52] was employed for the C, H, O and N atoms, while the Stuttgart RLC ECP basis set [53] was used for the Zn atom.

Molecular Dynamics Simulation
The crystallographic structure optimised at the QM/MM level and with the proper protonation (deprotonated GLU271, protonated GLU296) was chosen as the initial structure for the MD simulation. Protein parameters and charges were obtained from the ff14SB [37] force field from AMBER [35]. The substrate was parametrised using the parmchk2 module from AmberTools [54]. It was optimised at the B3LYP/6-31G(d) level, and then Merz-Kollman RESP charges [39] were calculated. Zn parameters were obtained following the Seminario Method [55] and the bonded model from AMBER using MCPB.py [56]. All residues within 3 Å around Zn were included in the model (HIS295, HIS299, GLU318, WAT1, WAT2). The model has not been optimised since the global structure corresponds to a QM/MM minimum, but frequencies were calculated as well as Merz-Kollman RESP charges [39] at the B3LYP/6-31G(d) level. MM parameters were derived from the frequencies' calculation.
AMBER's tLEAP module was used to generate an orthorhombic box of pre-equilibrated TIP3P [31] water molecules, including 9 Na + ions, to neutralise the protein's charge. Thus, a resulting box of 78.7 × 116.5 × 94.7 Å 3 containing almost 72,000 atoms, around 9700 of them coming from the protein, was obtained. Bonds between substrate's epoxide and WAT1 s H and substrate's epoxide and H of TYR383 were manually added by adding bond force constants of 25 kcal/mol·Å 2 in order to describe the detected H-bonds between these atoms.
MD calculations were performed in GPUs using the PMEMD module [57,58] in its CUDA version from the AMBER18 package.
Three initial MM minimisations were performed, combining the steepest descent method and the conjugate gradient method for 5000 steps each. In the first one, the protein, cofactors and substrate were kept frozen by applying a restraint of 150.0 kcal/mol·Å 2 . In the second one, only the backbone was restrained, while in the third minimisation, all the system was set free to move. A cut-off of 9 Å was applied in the three minimisations. Then, the system was heated from 0 to 300 K in steps of 30 K for 200 ps. SHAKE algorithm [59] was deactivated for the H of WAT and of TYR383, which are bonded to the epoxide. A restraint of 5 kcal/mol·Å 2 was applied to the backbone, and a cut-off of 9 Å was used.
Once heated, an NPT step of 5 ns was performed at 1 bar and at 300 K with a restraint of 5 kcal/mol·Å 2 applied to the backbone, a cut-off of 9 Å and with the SHAKE algorithm deactivated for the hydrogens bonded to the substrate's epoxide. A density of 1.0025 g/mL was achieved, and so the system was equilibrated. The temperature was controlled by Langevin dynamics [60] while the pressure was adjusted by the Berendsen barostat [61].
Finally, an NVT equilibration step was performed before the NVT production step. Both share the same configuration: 300 K, no restraints and no SHAKE algorithm for hydrogens bonded to substrate's epoxide. A total of 10 ns of equilibration were calculated followed by 150 ns of production, from which the last 100 ns have been selected for further analysis.

Conclusions
LTB 4 is a very potent lipid inflammatory mediator involved in acute and chronic inflammatory diseases. The chirality of the 12R-hydroxyl group and the 6Z,8E,10E configuration of the conjugated triene moiety are key features for its bioactivity. LTB 4 is obtained when LTA 4 is hydrolysed by the enzyme LTA 4 H. In this paper, our Molecular Dynamics simulations and Quantum Mechanics/Molecular Mechanics calculations allowed us to theoretically capture LTA 4 H in action, revealing the molecular details of the complete step-by-step mechanism of that enzyme reaction, in good agreement with the available experimental results.
The first step of the reaction consists of the opening of the epoxide of LTA 4 . In a second, well-separated step, a water molecule is added to C 12 of LTA 4 , while one of its hydrogen atoms is fully transferred to ASP375. The set of these two orthogonal steps constitutes a rather unusual 1,7-nucleophilic substitution through a clear S N 1 mechanism.
In LTA 4 , the bond C 6 -C 7 has a configuration rather similar to pro-cis, although intermediate between the pro-cis and the pro-trans ones. As a result of the epoxide opening, the bond C 6 -C 7 evolves either to a cis (Z) or to a trans (E) configuration. The triene moiety must twist to allow the bond C 6 -C 7 to adopt its cis (Z) or trans (E) configuration, in such a way that it exposes the R face of C 12 (by the side opposite to the Zn atom) or the S face of C 12 (by the same side as the Zn atom), respectively, to the water addition. Only the cis (Z) configuration will lead to LTB 4 (5S,12R-dihydroxy-6Z,8E,10E,14Z-eicosatetraenoic acid), whereas the trans (E) configuration would eventually form its stereoisomer 12S-6-trans-LTB 4 (5S,12S-dihydroxy-6E,8E,10E,14Z-eicosatetraenoic acid). Thus, the two main stereochemical features that are required for the bioactivity of LTB 4 appear to be closely linked.
In the gas phase, 12S-6-trans-LTB 4 is more stable than LTB 4 . Inside LTA 4 H, the epoxide opening (first step) to give the trans (E) configuration of the bond C 6 -C 7 turns out to be more favourable than the opening to the cis (Z) configuration, both kinetically and thermodynamically. However, from here on, the cis energy profile becomes clearly below the trans energy profile. This can be explained through noncovalent π-π stacking interactions between the triene moiety and TYR267 and, especially, between the triene moiety and TYR378. Both tyrosines wrap the triene system along the whole reaction.
Finally, in the third step, the final transfer of a proton from a water molecule to the oxygen atom bonded to C 6 (the former oxygen atom of the epoxide) takes place, thus forming LTB 4 . This proton transfer does not occur in the first step of the reaction to trigger an S N 1 acid-induced opening of the epoxide ring This is an excellent example of how the role of an enzyme is not only to accelerate the reaction rate, but to govern the stereochemistry of the product, making it possible for the bioactive product (LTB 4 in this case) to be formed. Funding: This study was supported by a grant (PID2020-113764GB-I00) from the Spanish "Ministerio de Ciencia e Innovación". We also acknowledge CSUC for computational facilities.