International Journal of Molecular Sciences Article Theoretical Characterization of the Step-by-Step Mechanism of Conversion of Leukotriene A4 to Leukotriene B4 Catalysed by the Enzyme Leukotriene A4 Hydrolase Miquel Canyelles-Niño 1,2 , Àngels González-Lafont 1,3 and José M. Lluch 1,3, * 1 Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain; miquel.canyelles@uab.cat (M.C.-N.); angels.gonzalez@uab.cat (À.G.-L.) 2 Biochemize SL, Carrer de Zamora, 45, 08005 Barcelona, Barcelona, Spain 3 Institut de Biotecnologia i de Biomedicina (IBB), Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain * Correspondence: josemaria.lluch@uab.cat Abstract: LTA4 H 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, LTA4 H 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 SN 1 mechanism, the epoxide Citation: Canyelles-Niño, M.; opens and the triene moiety of the substrate twists in such a way that the bond C6 -C7 adopts its cis (Z) González-Lafont, À.; Lluch, J.M. configuration, thus exposing the R face of C12 to the addition of a water molecule hydrogen-bonded Theoretical Characterization of the to ASP375. Thus, the two stereochemical features that are required for the bioactivity of LTB4 appear Step-by-Step Mechanism of Conversion of Leukotriene A4 to to be closely related. The noncovalent π-π stacking interactions between the triene moiety and two Leukotriene B4 Catalysed by the tyrosines (TYR267 and, especially, TYR378) that wrap the triene system along the whole reaction Enzyme Leukotriene A4 Hydrolase. explain the preference for the cis configuration inside LTA4 H. Int. J. Mol. Sci. 2022, 23, 3140. https://doi.org/10.3390/ Keywords: leukotriens; leukotriene A4 hydrolase; enzyme catalysis; QM/MM calculations; molecular ijms23063140 dynamics simulations; proinflammatory lipid mediators Academic Editor: Francesco Mallamace Received: 17 February 2022 1. Introduction Accepted: 11 March 2022 Nowadays, it is recognised that inflammatory-based human diseases represent the Published: 15 March 2022 leading causes of present-day morbidity and mortality worldwide [1]. Infection of the host, Publisher’s Note: MDPI stays neutral tissue injury or surgical trauma trigger the release of proinflammatory lipid and protein with regard to jurisdictional claims in mediators by cells in the site of challenge. These chemical mediators act as chemoattractants published maps and institutional affil- to recruit neutrophiles. This initial inflammatory response occurs to protect the host and iations. 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 Copyright: © 2022 by the authors. does not work well, the acute inflammation evolves to chronic inflammation [3–7]. Chronic Licensee MDPI, Basel, Switzerland. inflammatory diseases are the most significant cause of death in the world today [1]. This article is an open access article The lipid mediator Leukotriene B4 (LTB4 , 5S,12R-dihydroxy-6Z,8E,10E,14Z-eicosa- distributed under the terms and conditions of the Creative Commons tetraenoic acid) (Figure 1a) exerts one of the most potent chemotactic effects on polymor- Attribution (CC BY) license (https:// phonuclear leukocytes in the onset phase. In addition, excessive formation of LTB4 can creativecommons.org/licenses/by/ be related to the maintenance of chronic inflammatory diseases [8]. The biosynthesis of 4.0/). LTB4 from arachidonic acid (AA) is the result of the action of several proteins. To start, Int. J. Mol. Sci. 2022, 23, 3140. https://doi.org/10.3390/ijms23063140 https://www.mdpi.com/journal/ijms from arachidonic acid (AA) is the result of the action of several proteins. To start, the 5- lipoxygenase-activating protein (FLAP), an integral membrane protein, facilitates the transfer of the substrate AA to the enzyme 5-lipoxygenase (5-LOX) [9–11]. 5-LOX pro- duces 5S-hydroperoxy-6E,8Z,11Z,14Z-eicosatetraenoic acid (5-HpETE) as a result of the hydroperoxidation Int. J. Mol. Sci. 2022, 23, 3140 of C5 of AA. Then, 5-LOX catalyses the formation of the highly labile 2 of 23 epoxide intermediate Leukotriene A4 (LTA4, 5S,6S-epoxy-7E,9E,11Z,14Z-eicosate-traenoic acid) (Figure 1b) by dehydration of 5-HpETE. FLAP regulates the 5-LOX activity, increas- ing the efficiencythe of 5-lipoxygenase-activating the formation of LTA4 fromprotein5-HpETE [8,12–14]. (FLAP), Finally,membrane an integral after havingprotein, facilitates been released fromthe 5-LOX, transferLTA 4 is hydrolased into LTB4 by the enzyme leukotriene A4 of the substrate AA to the enzyme 5-lipoxygenase (5-LOX) [9–11]. 5-LOX hydrolase (LTA4produces H) (EC 3.3.2.6) [8]. 5S-hydroperoxy-6E,8Z,11Z,14Z-eicosatetraenoic acid (5-HpETE) as a result of the LTA4H is a hydroperoxidation monomeric bifunctional of C5 zinc metalloenzyme. of AA. It is detected Then, 5-LOX catalyses in almost all the formation of the highly labile mammalian cells,epoxide and it has both aminopeptidase activity and epoxide hydrolase activity. intermediate Leukotriene A4 (LTA4 , 5S,6S-epoxy-7E,9E,11Z,14Z-eicosate-traenoic In this paper, weacid) will focus on1b) (Figure thebysecond one. This dehydration epoxide hydrolase of 5-HpETE. reaction FLAP regulates is 5-LOX the interest- activity, increasing ing because involves the stereoselective introduction of one water molecule to the carbon the efficiency of the formation of LTA4 from 5-HpETE [8,12–14]. Finally, after having been backbone (C12) ofreleased LTA4 infroma position 5-LOX,several LTA4 ismethylene hydrolasedunits intoaway LTB4 from by thethe epoxyleukotriene enzyme moiety A4 hydrolase (carbons C5 to C6(LTA ) [15].4 H) (EC 3.3.2.6) [8]. (a) (b) Figure Figure 1. Structures 1. Structures corresponding corresponding to LTB to LTB (a) and LTA4 (b). 4 (a) and LTA4 (b).4 LTA4 H is a monomeric bifunctional zinc metalloenzyme. It is detected in almost all A high-resolution crystal structure of human LTA4H in the complex with the com- mammalian cells, and it has both aminopeptidase activity and epoxide hydrolase activity. In petitive inhibitor bestatin revealed that the protein includes three domains, an N-terminal this paper, we will focus on the second one. This epoxide hydrolase reaction is interesting domain, a Zn-containing catalytic domain and an α-helical C-terminal domain. The active because involves the stereoselective introduction of one water molecule to the carbon site is placed in a deep cleft in between the domains, where the Zn2+ is coordinated to backbone (C12 ) of LTA4 in a position several methylene units away from the epoxy moiety HIS295, HIS299 and one carboxylic oxygen of GLU318 [16,17]. Because of the LTA4 insta- (carbons C5 to C6 ) [15]. bility (half-life of about 10 s at neutral pH) high-resolution crystal structures of LTA4H A high-resolution crystal structure of human LTA4 H in the complex with the compet- complexed to LTA4 were not determined until recently by Haeggström and coworkers itive inhibitor bestatin revealed that the protein includes three domains, an N-terminal [18]. They described domain,six different structures a Zn-containing of human catalytic domain LTA 4H from five distinct crystal and an α-helical C-terminal domain. The active forms, finding two site is placed in a deep cleft in between theseveral conformational states of the enzyme and conformations domains, where the Zn of 2+ theis coordinated to substrate LTA4, and showing that LTA 4H undergoes domain movements. From this struc- HIS295, HIS299 and one carboxylic oxygen of GLU318 [16,17]. Because of the LTA4 insta- tural informationbility and previous (half-life site-directed of about 10 smutagenesis at neutral pH) experiments [15,19,20], high-resolution Haegg- crystal structures of LTA4 H ström and coworkers [18] propose complexed a mechanism to LTA4 were for theuntil not determined epoxide hydrolase recently reaction. and by Haeggström In coworkers [18]. short, the C1 carboxylate group of LTA 4 anchors to ARG563, while Zn2+ and TYR383 are They described six different structures of human LTA4 H from five distinct crystal forms, coordinated to the epoxide finding two oxygen. A catalyticstates conformational waterofmolecule the enzyme(WAT1), polarised and several by Zn2+ conformations of the sub- and GLU271, transfers a proton strate LTA to produce an SN1 acid-induced opening of the epoxide 4 , and showing that LTA4 H undergoes domain movements. From this structural ring. The positiveinformation charge of the andresulting previous carbonium ion mutagenesis site-directed is delocalisedexperiments over the triene sys- [15,19,20], Haeggström tem, with the C6-C and7 bond in a pro-cis coworkers configuration [18] propose stabilised a mechanism by the for the L-shape epoxide of the tight hydrolase reaction. In short, the pocket. Finally, aCsecond water molecule, activated by ASP375, adds to C12 to 1 carboxylate group of LTA4 anchors to ARG563, while Zn 2+ generate and TYR383 an are coordinated to the epoxide oxygen. A catalytic water molecule (WAT1), polarised by Zn2+ and GLU271, transfers a proton to produce an SN 1 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 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 LTB4 [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]. Int. J. Mol. Sci. 2022, 23, 3140 3 of 23 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–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 LTB4 . In this paper, we combine Molecular Dynamics (MD) simulations and Quantum Me- chanics/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 A4 to leukotriene B4 catalysed by the enzyme leukotriene A4 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 LTB4 . 2. Results and Discussion 2.1. 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 LTA4 :LTA4 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 LTA4 :LTA4 H Michaelis complex is pictured in Figure 3. A scheme of the main noncovalent interactions between the substrate LTA4 and the enzyme LTA4 H in this optimised Michaelis complex is shown in Figure 4. LTA4 fits into the L-shaped hydrophobic cavity of the LTA4 H. 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 LTA4 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 LTA4 , GLU271 (hydrogen-bonded to WAT1 and WAT2) and GLU296 (linked to WAT1 through a hydrogen bond). Int. Int. J. J. Mol. Mol. Sci. 2022, 23, Sci. 2022, 23, 3140 x FOR PEER REVIEW 44 of of 23 27 Figure 2. Figure 2. QM/MM QM/MM partition. partition. QM QM atoms atoms are are depicted depicted in in red. red. The The boundary boundary between the QM between the QM and and MM regions is indicated by blue wavy lines. MM regions is indicated by blue wavy lines. 2.2. Molecular Dynamics Simulation In the above section we have analysed the structure of the optimised LTA4 :LTA4 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 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 Å Int. J. Mol. Sci. 2022, 23, 3140 5 of 23 Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 5 of 27 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. Figure 3. Complete view of the QM/MM-optimised LTA4 :LTA4 H Michaelis complex (the solvation Figure 3. Complete waters view of are not shown forthe theQM/MM-optimised sake of clarity). LTALTA 4:LTA 4H has 4H Michaelis three domainscomplex (the that are solvation pictured here waters are not colours: in different shown forThethe sake of clarity). N-terminal domain LTA (in4H has three green), domainsdomain the catalytic that are(in pictured yellow)here andinthe different colours: C-terminal The (in domain N-terminal domain violet). The (in green), substrate LTA4 isthe catalytic using represented domain (in yellow) a stick model. and the C- terminal domain (in violet). The substrate LTA4 is represented using a stick model. 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 |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 LTB4, as will be explained below. 2.3. 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 pos- sible. To unravel this point, a two-dimensional potential energy surface was built (see Figure 9) using a reaction coordinate rc1 for the epoxide ring opening (rc1 = d(C6(LTA4) − Oepox(LTA4)), that is, the epoxide breaking bond length, and a reaction coordinate rc2 for the water addition to C12 (rc2 = (d(O(WAT3) − H(WAT3)) − d(O(WAT3) − C12(LTA4))), that is, the difference be- tween 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). Int. J.J. Mol. Int. Mol. Sci. Sci. 2022, 2022, 23, 23, 3140 x FOR PEER REVIEW 66 of 27 of 23 Figure 4. Figure 4. Scheme Schemeofofthe main the noncovalent main interactions noncovalent (in dashed interactions lines) (in dashed between lines) the substrate between LTA4 the substrate (in black) and the enzyme LTA4H (in green) in the optimised Michaelis complex and the Zn envi- LTA4 (in black) and the enzyme LTA4 H (in green) in the optimised Michaelis complex and the Zn ronment. The dashed orange lines indicate π-π stacking interactions. environment. The dashed orange lines indicate π-π stacking interactions. LTA As 4 fits into the L-shaped hydrophobic cavity of the LTA4H. The substrate docks to shown in Figure 9, the reaction takes place in two well-separated steps. Firstly, the the protein C6 − Oepox bond through the interaction progressively between breaks, the polar thus opening thehead of thering. epoxide LTAIn 4 and ARG563 and the intermediate INT1, the epoxide is already open (the distance C6 − Oepox is 2.3 Å), but thethe LYS565 from the protein, which are located in the depth of the cavity. Thus, substrate value of the is placed in a head-first2 conformation. Moreover, this fitting conformation reaction coordinate rc has not significantly changed yet. From here on, in the second lets the epoxide step oxygen the moveC6near distance the Znremains − Oepox environment, while it places quite invariant, whileCa12clear near evolution ASP375 (with of theWAT3 be- reaction tween them). coordinate rc2 The epoxide occurs, is not only indicating stabilised the water by being (WAT3) in the addition to second C12 andcoordination the breakagesphere of an of the Zn atom, but also thanks to two hydrogen bonds O-H bond in this water molecule and finally reaching the intermediate INT2. from TYR383 and WAT1, the cat- Thus, the alytic water which 1,7-nucleophilic is also coordinated substitution to Zn. clearly takes With place this conformation, through TYR267 and an SN 1 mechanism. As aTYR378 matter arefact, of the ittwo is an residues extremeclosest case oftoSNthe triene moiety 1 reaction, becauseof the substrate. The the minimum energymain interactions reaction paths between the substrate LTA corresponding to each step appear to be roughly orthogonal (each one parallelcomplex 4 and the enzyme LTA 4H in this optimised Michaelis to one are schematised coordinate axis) in inthe Figure 4. The role ofpotential two-dimensional all these energy residues is to limit surface shown theinavailable Figure 9.space Any and give aofshape deviation to the cavity. the mechanism TYR267 towards theand TYR378 region neararethekey residues diagonal of the of that selectivity energy surfaceof the reaction too (see below). (that is, any approach to the SN 2 mechanism) implies the penetration in high-energy The Zn ion is the central atom of the catalytic domain of the protein. In good agree- ment with the crystallographic structures [18], in the optimised structure it is pentacoor- dinated to two HISs (HIS295, HIS299), GLU318, WAT1 and WAT2 with an octahedral ge- Int. J. Mol. Sci. 2022, 23, 3140 ometry (distance RMSD compared to the ideal octahedron: 0.490 Å) [28], so a 7vacancy of 23 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 ini- regions, which turns out to be forbidden. An NBO charge analysis in INT1 gives a charge tially contains the epoxide of LTA4, GLU271 (hydrogen-bonded to WAT1 and WAT2) and of +0.71 a.u. delocalised over the triene system (from C6 to C12 ) of the substrate, which GLU296 (linked to WAT1 through a hydrogen bond). confirms the existence of the carbocation required in the intermediate of an SN 1 mechanism. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 8 of 27 Figure 5. Zn environment in the optimised structure. Figure 5. Zn environment in the optimised structure. 2.2. 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 sys- tem. To unravel this point, we analysed 100 ns extracted from a Molecular Dynamics sim- ulation 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 Distances all the snapshots ASP375-WAT (blue (96%Cof-WATthe generated snapshots) theinclude a bridg- Figure Figure 6. 6. Distances ASP375-WAT (bluepoints) points)and and 12 C12-WAT (orange points) (orange along points) along Molecular the Molecular ing water Dynamics molecule closer simulationfor forthe than the LTA 3 LTA44:LTA Å and 5 Å to ASP375 and C 12, respectively. Even as many Dynamics simulation :LTA44HHMichaelis Michaeliscomplex. complex.WATWATstands for afor stands water molecule a water molecule as 58% placed placed of theASP375 between between generatedand snapshots ASP375and CC12 12.. 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 mole- cule, Onactivated the otherby ASP375, hand, appears to the evolution of be ready to distances attackinCFigure shown 12, thus forming a hydroxyl 7 confirms that the group. epoxide of LTA4 maintains two hydrogen bonds with TYR383 and WAT1, and that it re- mains in the second coordination sphere of Zn ion along the Molecular Dynamics simu- lation. 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 re- Int. J. Mol. Sci. 2022, 23, 3140 mains in the second coordination sphere of Zn ion along the Molecular Dynamics 8 of 23simu- lation. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 9 of 27 Figure Figure 7. Distancesepoxide-WAT1 7. Distances epoxide-WAT1 (blue (blueline), line),epoxide-TYR383 epoxide-TYR383 (orange line)line) (orange and and epoxide-Zn (green(green epoxide-Zn line) along the Molecular Dynamics simulation for the line) along the Molecular Dynamics simulation for the LTA LTA :LTA 4 4:LTA 4 H Michaelis complex. 4H Michaelis complex. 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 simula- tion. 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 pro- trans configuration |𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 𝑎𝑛𝑔𝑙𝑒 𝐻6 𝐶6 𝐶7 𝐻7| 90°. The evolution of this di- hedral angle along the complete enzyme reaction is a key point for the formation of LTB4, as will be explained below. Figure 8. Dihedral angle H -C -C -H along the Molecular Dynamics simulation for the LTA :LTA H Figure 8. Dihedral angle H66-C66-C77-H77along the Molecular Dynamics simulation for the4LTA44:LTA 4H Michaelis complex. The pro-cis configuration (orange line) corresponds to |dihedral angle| ≤ 90◦ , Michaelis complex. The pro-cis configuration (orange line) corresponds to |𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 𝑎𝑛𝑔𝑙𝑒| 90°, whereas |dihedral angle| ≥ 90◦ indicates a pro-trans configuration (blue line). whereas |𝑑𝑖ℎ𝑒𝑑𝑟𝑎𝑙 𝑎𝑛𝑔𝑙𝑒| 90° indicates a pro-trans configuration (blue line). 2.3. 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 open- ing ( 𝑟 = d(C6(LTA4)−Oepox(LTA4)), that is, the epoxide breaking bond length, and a reac- Int. Int. J.J. Mol. Mol. Sci. Sci. 2022, 2022, 23, 23, x3140 FOR PEER REVIEW 10 9of of 27 23 Figure Figure 9. Two-dimensional Two-dimensional potential potential energy energy surface surface as as a function function of of the epoxide epoxide ring ring opening opening and and water water addition reaction coordinates defined in the text. Isoenergetic lines correspond to intervals addition reaction coordinates defined in the text. Isoenergetic lines correspond to intervals of of 55 kcal/mol. kcal/mol. As shown To our in Figure surprise, 9, the reaction a structural analysis takes of theplace in two well-separated intermediate INT2 showssteps. that the Firstly, con- the C − O bond progressively breaks, thus opening the jugated triene moiety has a 6E,8E,10E configuration that would not lead to LTB4 but to 6 epox epoxide ring. In the intermedi- ate INT1, the epoxide its stereoisomer is already4open 12S-6-trans-LTB (the distance C6 − Oepox is 2.3 Å), but the value (5S,12S-dihydroxy-6E,8E,10E,14Z-eicosatetraenoic of acid). the reactionthe Following evolution𝑟of the coordinate has structures not significantly along the changed yet. From here SN 1 mechanism shown on,ininFigure the second 9, we step found the distance that C6 − Oepox an unexpected remains rotation quite of the C6 -Cinvariant, 7 bond from while theainitial clear evolution rather pro-cis of the reac- configu- tion ration in LTA4 to𝑟a pro-trans coordinate occurs, indicating configuration the water in INT1 (WAT3) occurs,addition which leads to Cto 12 and the 6E theinadequate breakage of an O-H bond configuration of in C6this -C7 water in INT2. molecule To analyseand this finally fact,reaching we builtthe intermediate INT2. a two-dimensional Thus, potential energy surface (see substitution Figure 10) asclearlya function 3 the 1,7-nucleophilic takesof the through place rotation an of the SN1 C 6 -C7 bond (r mechanism. Asc = a dihe- mat- dralofangle ter fact, H -Can it 6is 6 -Cextreme 7 -H7 ) and theof case epoxide ring opening SN1 reaction, because (rc1the = d(C − Oepoxreaction 6 (LTA4 ) energy minimum (LTA4 )). Starting paths from the reactant corresponding to each R, step the Dijsktra appear to algorithm be roughly [29]orthogonal localises two (eachminimum one parallel energy to reaction one paths (grey coordinate axis) in andthewhite points) on the two-dimensional surface energy potential shown surface in Figure 10, and shown the mini- in Figure 9. mum (R, INT1trans and INT1cis) and maximum (TS1trans and Any deviation of the mechanism towards the region near the diagonal of that energy sur- TS1cis) energy points along them. face One (that is,path any involves approachthe pro-trans to the (grey points) SN2 mechanism) configuration implies the penetration(as a matter of fact, it in high-energy corresponds to the first step regions, which turns out to be forbidden. of the S 1 mechanism shown in Figure N An NBO charge analysis in INT1 gives a charge 9). The other one maintains of +0.71 a.u. a pro-cis delocalisedconfiguration (whitesystem over the triene points)(fromto reach C6 to anCintermediate 12) of the substrate, INT1ciswhich with the adequate confirms pro-cis structure. the existence The two paths of the carbocation requiredare competitive, in the intermediate with the of pro-trans an SN1 mecha- being the one that implies the lowest potential energy barrier. This explains why only the trans nism. configuration To our surprise, was obtained in the a structural two-dimensional analysis potential of the intermediate energy INT2 shows surface that theshown in conju- Figure 9, where the rotation of the gated triene moiety has a 6E,8E,10E configuration C -C 6 7 bond was not included to define that would not lead to LTB4 but to its the reaction coordinate. This stereoisomer result also agrees 12S-6-trans-LTB with our B3LYP/6-31G(d) calculation that shows that in 4 (5S,12S-dihydroxy-6E,8E,10E,14Z-eicosatetraenoic acid). the gas phase, Following 12S-6-trans-LTB the evolution is 0.3 kcal/mol of the4structures alongmore the SNstable than LTBshown 1 mechanism 4 in terms of potential in Figure 9, we energy. However, the experimental result is that LTA 4 H converts found that an unexpected rotation of the C6-C7 bond from the initial rather pro-cis config- LTA 4 into LTB 4 , that is, with a 6Z,8E,10E configuration for the conjugated triene moiety. uration in LTA4 to a pro-trans configuration in INT1 occurs, which leads to the 6E inade- quateAn especially interesting configuration of C6-C7 in point INT2. is Tothat the configuration analyse this fact, we built of the C6 -C7 bond in INT1 a two-dimensional po- determines the stereochemistry of the ingoing water tential energy surface (see Figure 10) as a function of the rotation of the into C 12 . In Figure C6-C11, we have 7 bond (𝑟 = displayed the position of WAT3 with respect to the plane generated by the triene moiety dihedral angle H6-C6-C7-H7) and the epoxide ring opening (𝑟 = d(C6(LTA4)−Oepox(LTA4)). (from C to C ) in the reactant LTA4 (pro-cis, yellow disk) in INT1trans (pro-trans, brown Starting6from12the reactant R, the Dijsktra algorithm [29] localises two minimum energy disk) and in INT1cis (pro-cis, violet disk). It can be clearly seen that the face of attack reaction paths (grey and white points) on the surface shown in Figure 10, and the mini- of WAT3 in the pro-trans configuration is just the opposite than in the case of the pro-cis mum (R, INT1trans and INT1cis) and maximum (TS1trans and TS1cis) energy points configuration in INT1, leading to the 12S-hydroxyl or the 12R-hydroxyl configurations, 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 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 con- figuration 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. Int. J. Mol. Sci. 2022, 23, 3140 10 of 23 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 respectively. Thus, the two 6Z,8E,10E configuration formain stereochemical the conjugated features triene that are required for the bioactivity moiety. of LTB4 appear to be closely linked. Figure 10. Figure 10. Two-dimensional potential energy surface as a function function of of the the rotation rotationof ofthe theCC6-C 6 -C 77bond bond Int. J. Mol. Sci. 2022, 23, x FOR PEERand REVIEW and the the epoxide epoxide ring-opening reaction reaction coordinates coordinates defined definedin inthe thetext. text.Isoenergetic Isoenergeticlines 12 of 27 linescorrespond correspond tointervals to intervals of of 55 kcal/mol. kcal/mol. 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 dis- played 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 bioac- tivity of LTB4 appear to be closely linked. Figure 11. Position of WAT3 with respect to the plane generated by the triene moiety in the reactant 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). 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 At this point, it is clear that a complete and comparative study of the pro-cis and pro-trans reaction paths is needed to understand the mechanism of this enzyme reaction. the pro-trans reaction paths is needed to understand the mechanism of this enzyme re- To this aim, firstly we fully optimised the structures INT1cis, INT1trans, TS1cis and action. To this aim, firstly we fully optimised the structures INT1cis, INT1trans, TS1cis TS1trans, which appear in Figure 10, in order to locate the corresponding stationary struc- and TS1trans, which appear in Figure 10, in order to locate the corresponding stationary tures (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 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 Int. J. Mol. Sci. 2022, 23, 3140 pro-trans reaction paths is needed to understand the mechanism of this enzyme reaction. 11 of 23 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 struc- tures (minima structures (minima and and transition-state structures) transition-state in the structures) incomplete potential the complete hypersurface potential hypersur- of the reaction. Then, for each reaction path, we followed the reaction face of the reaction. Then, for each reaction path, we followed the reaction coordinates coordinates d(O(WAT3)-H(WAT3))-d(O(WAT3)-C12(LTA d(O(WAT3)-H(WAT3))-d(O(WAT3)-C (LTA)) 4)) and d(O(WAT1)-H(WAT1))- 12 4 and d(O(WAT1)-H(WAT1))-d(H(WAT1)- d(H(WAT1)-O epox) to describe, respectively, the water addition (WAT3) to C12 and the final Oepox ) to describe, respectively, the water addition (WAT3) to C12 and the final protonation ofprotonation the oxygen of atomthe bonded oxygen toatom bonded C6 (the to Coxygen former 6 (the former oxygen atom of the epoxide) by atom of the epoxide) by WAT1 and to WAT1the locate and to locate the corresponding corresponding stationary stationary structures. Thestructures. The relative relative potential potential energies for allener- the gies for all the stationary points we localised stationary points we localised are shown in Figure 12. are shown in Figure 12. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 13 of 27 Figure 12. Diagram of potential energies of the stationary structures of the complete mechanism of Figure conversion of LTA4 to LTB412. Diagram (pro-cis of potential reaction path,energies in blue) ofcatalysed the stationary structures by LTA 4 H andof the complete of the mechanism of competitive conversion of LTA4 to LTB4 (pro-cis reaction path, in blue) catalysed by LTA4H and of the competi- pro-trans reaction path (in red). tive pro-trans reaction path (in red). We focused firstlyWe onfocused the pro-cis reaction path. In Figure 13 we show a succession of firstly on the pro-cis reaction path. In Figure 13 we show a succession of pictures of the stationary structures along this path, pictures of the stationary structures alongintending this path, to capture intending toLTA 4 H LTA capture in action. 4H in action. For the sake of clarity, For theinsake Figures S1–S3 of clarity, we display in Figures S1–S3 representations of thoseofstationary we display representations those stationary structures focusedstructures on the region focusedwhere on thethe corresponding region step occurs where the corresponding in occurs step each case. in each case. R Figure 13. Cont. Int. J. Mol. Sci. 2022, 23, 3140 12 of 23 R Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 14 of 27 TS1cis INT1cis TS2cis Figure 13. Cont. Int. J. Mol. Sci. 2022, 23, 3140 13 of 23 Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 15 of 27 INT2cis Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 16 of 27 TS3cis Pcis Figure Figure 13. Stationary 13. Stationary structures of thestructures completeofmechanism the complete of mechanism conversionof conversion of LTA4 toof LTA LTB44to(pro-cis LTB4 (pro-cis reaction path) catalysed by LTA4H: R, TS1cis, INT1cis, TS2cis, INT2cis, TS3cis and Pcis. reaction path) catalysed by LTA4 H: R, TS1cis, INT1cis, TS2cis, INT2cis, TS3cis and Pcis. 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 Int. J. Mol. Sci. 2022, 23, 3140 14 of 23 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 moves away from INT1cis. Conversely, the Zn-Oepox 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 Oepox -C5 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 SN 1 mechanism. Very interestingly, the dihedral angle H6 -C6 -C7 -H7 goes from −65.6◦ in R (a value close to an intermediate situation of the bond C6 -C7 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 C12 (by the side opposite to the Zn atom, see Figure S2) at INT1cis. Then, WAT3 adds to C12 (the distance O(WAT3)-C12 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 C6 (the former oxygen atom of the epoxide, Oepox ) takes place, while the distance Zn-Oepox clearly increases through a potential energy barrier of 14.4 kcal/mol, thus forming LTB4 . Note that Haeggström and coworkers [18] have proposed that the catalytic water WAT1 transfers a proton to promote an SN 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 C6 -C7 bond for the stationary structures of the complete mechanism of conversion of LTA4 to LTB4 (pro-cis reaction path) catalysed by LTA4 H. Pro-cis Reaction Path R TS1 INT1 TS2 INT2 TS3 P Zn-Oepox 3.39 2.94 1.99 1.97 1.95 2.45 2.96 Oepox -C5 1.45 1.41 1.39 1.39 1.39 1.41 1.43 Oepox -C6 1.47 1.92 2.40 2.41 2.42 2.37 2.39 H(WAT1)-Oepox 1.75 1.65 2.42 2.46 2.49 1.39 1.02 O(WAT1)-H(WAT1) 0.99 1.00 0.97 0.97 0.98 1.09 1.59 O(WAT3)-C12 3.84 3.69 3.69 1.82 1.53 1.53 1.53 H(WAT3)-O(ASP375) 1.71 1.69 1.67 1.36 1.03 1.03 1.03 O(WAT3)-H(WAT3) 0.99 0.99 0.99 1.11 1.53 1.53 1.53 H6 -C6 -C7 -H7 −65.6 −15.9 −7.1 1.5 1.5 −0.5 −0.7 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. H(WAT3)-O(ASP375) 1.71 1.69 1.67 1.36 1.03 1.03 1.03 O(WAT3)-H(WAT3) 0.99 0.99 0.99 1.11 1.53 1.53 1.53 H6-C6-C7-H7 −65.6 −15.9 −7.1 1.5 1.5 −0.5 −0.7 Let us describe the pro-trans reaction path now. In Figure 14, we show a succession Int. J. Mol. Sci. 2022, 23, 3140 of pictures of the stationary structures along this path. For the sake of clarity,15inofFigures 23 S4–S6, we display representations of those stationary structures focused on the region where the corresponding step occurs in each case. Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 18 of 27 TS1trans INT1trans TS2trans Figure 14. Cont. Int. J. Mol. Sci. 2022, 23, 3140 16 of 23 Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 19 of 27 INT2trans Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 20 of 27 TS3trans Ptrans Figure Figure 14. Stationary 14. Stationary structures of thestructures complete of mechanism the complete mechanism of the pro-trans of the pro-trans reactionreaction path catalysed path catalysed by LTA4H: TS1trans, INT1trans, TS2trans, INT2trans, TS3trans and Ptrans. by LTA4 H: TS1trans, INT1trans, TS2trans, INT2trans, TS3trans and Ptrans. 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- Int. J. Mol. Sci. 2022, 23, 3140 17 of 23 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-Oepox bond is significantly shortened from 3.39 Å in R to 1.98 Å at INT1trans. The epoxide is fully opened at INT1trans of this SN 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 H6 -C6 - C7 -H7 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 C12 (by the same side as the Zn atom, see Figure S5) at INT1trans. Then, WAT3 adds to C12 (the distance O(WAT3)-C12 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 C6 takes place, while the distance Zn-Oepox clearly increases through a potential energy barrier of 17.8 kcal/mol, thus forming 12S-6-trans-LTB4 , the stereoisomer of LTB4 . Table 2. Most relevant distances (in Å) and dihedral angle (in degrees) defining the E/Z character of the C6 -C7 bond for the stationary structures of the pro-trans reaction path catalysed by LTA4 H. Pro-trans Reaction Path R TS1 INT1 TS2 INT2 TS3 P Zn-Oepox 3.39 2.64 1.98 1.93 1.92 3.26 2.96 Oepox -C5 1.45 1.40 1.39 1.39 1.39 1.44 1.43 Oepox -C6 1.47 2.04 2.20 2.30 2.32 2.39 2.34 H(WAT1)-Oepox 1.75 1.66 3.00 3.12 3.12 0.97 1.02 O(WAT1)-H(WAT1) 0.99 1.00 0.98 0.98 0.98 2.39 1.59 O(WAT3)-C12 3.84 3.68 3.66 1.79 1.51 1.50 1.50 H(WAT3)-O(ASP375) 1.71 1.66 1.66 1.32 1.01 1.01 1.01 O(WAT3)-H(WAT3) 0.99 1.00 1.00 1.12 1.60 1.64 1.62 H6 -C6 -C7 -H7 294.4 196.9 183.7 178.1 177.5 180.7 181.5 A final point to discuss is why LTA4 H converts LTA4 into LTB4 instead of the stereoisomer 12S-6-trans-LTB4 , 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 LTA4 H formation of LTB4 is more favourable than the formation of 12S-6-trans-LTB4 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◦ ) Int. J. Mol. Sci. 2022, 23, 3140 18 of 23 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 C6 -C7 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-LTB4 and leading to LTB4 as a product of the reaction catalysed by LTA4 H. Indeed, those two tyrosines do not exist in the gas phase, where the 12S-6-trans-LTB4 is more stable than LTB4 . Our theoretical result agrees Int. J. Mol. Sci. 2022, 23, x FOR PEER REVIEW 22 with of 27 experimental data [21,22] showing that mutations of TY378 lead not only to LTB4 but also to products with a trans stereochemistry in the C6 -C7 bond. Anglesbetween 15. Angles Figure 15. betweenthe theplane plane ofofthethe triene triene moiety moiety of the of the substrate substrate andand the planes the planes of theofaro- the matic rings aromatic of TYR267 rings of TYR267andandTYR378 TYR378 forfor the thedifferent differentstationary stationarystructures structureslocated locatedalong along the the pro-cis pro-cis and and the the pro-trans pro-trans reaction reaction paths. paths. The The closest closest plane plane to to the the whole whole number number ofof C C and and H H atoms atoms of of the the triene triene system system was was taken taken as the plane as the of this plane of this triene triene moiety. moiety. 3. Materials Materials and Methods 3.1. Protein 3.1. Protein Setup Setup Numerous crystallographic Numerous crystallographic structures structures [18] [18] have have been been reported reported forfor human human LTALTA44H, H, some of them complexing a wide variety of inhibitors and a few of them containing some of them complexing a wide variety of inhibitors and a few of them containing the the LTA44 substrate LTA substrate but but including one or including one or more more mutations. Here, aa single-mutated mutations. Here, version of single-mutated version of human LTA 4 H (ASP375ASN) complexed with LTA 4 was used (PDB code: human LTA4H (ASP375ASN) complexed with LTA4 was used (PDB code: 5NI6) [18]. The 5NI6) [18]. The mutation was mutation was reverted reverted toto the the WT WT enzyme. The protein enzyme. The protein was was protonated protonated at at pH pH 7.0 7.0 using using PropKa3.0 [30] through a web interface (www.playmolecule.org, accessed on PropKa3.0 [30] through a web interface (www.playmolecule.org, access on 1 July 2019). 1 July 2019). Two versions of the protonation state were prepared: (GLU271, protonated GLU296) and Two versions of the protonation state were prepared: (GLU271, protonated GLU296) and (protonated GLU271, GLU296), leaving all the other protonation states untouched. A box (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 of TIP3P [31] water molecules was built around the protein for solvation. The box was up considering 10 Å from the outermost protein’s atom in each direction of the space. built up considering 10 Å from the outermost protein’s atom in each direction of the space. 3.2. 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- Int. J. Mol. Sci. 2022, 23, 3140 19 of 23 3.2. 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 (LTA4 ), 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 C4 to C13 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 C3 -C4 and C13 -C14 ). 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. 3.3. 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 Int. J. Mol. Sci. 2022, 23, 3140 20 of 23 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 frequen- cies’ 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 WAT10 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. 4. Conclusions LTB4 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 configu- ration of the conjugated triene moiety are key features for its bioactivity. LTB4 is obtained when LTA4 is hydrolysed by the enzyme LTA4 H. In this paper, our Molecular Dynamics simulations and Quantum Mechanics/Molecular Mechanics calculations allowed us to theoretically capture LTA4 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 LTA4 . In a second, well-separated step, a water molecule is added to C12 of LTA4 , 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 SN 1 mechanism. In LTA4 , the bond C6 -C7 has a configuration rather similar to pro-cis, although inter- mediate between the pro-cis and the pro-trans ones. As a result of the epoxide opening, the bond C6 -C7 evolves either to a cis (Z) or to a trans (E) configuration. The triene moiety must twist to allow the bond C6 -C7 to adopt its cis (Z) or trans (E) configuration, in such a way that it exposes the R face of C12 (by the side opposite to the Zn atom) or the S face of C12 (by the same side as the Zn atom), respectively, to the water addition. Only the cis (Z) configuration will lead to LTB4 (5S,12R-dihydroxy-6Z,8E,10E,14Z-eicosatetraenoic acid), Int. J. Mol. Sci. 2022, 23, 3140 21 of 23 whereas the trans (E) configuration would eventually form its stereoisomer 12S-6-trans-LTB4 (5S,12S-dihydroxy-6E,8E,10E,14Z-eicosatetraenoic acid). Thus, the two main stereochemical features that are required for the bioactivity of LTB4 appear to be closely linked. In the gas phase, 12S-6-trans-LTB4 is more stable than LTB4 . Inside LTA4 H, the epoxide opening (first step) to give the trans (E) configuration of the bond C6 -C7 turns out to be more favourable than the opening to the cis (Z) configuration, both kinetically and ther- modynamically. 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 C6 (the former oxygen atom of the epoxide) takes place, thus forming LTB4 . This proton transfer does not occur in the first step of the reaction to trigger an SN 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 (LTB4 in this case) to be formed. Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms23063140/s1. Author Contributions: Conceptualization, À.G.-L. and J.M.L.; investigation, M.C.-N.; writing— original draft preparation, M.C.-N., À.G.-L. and J.M.L.; writing—review and editing, À.G.-L. and J.M.L.; funding acquisition, À.G.-L. and J.M.L. All authors have read and agreed to the published version of the manuscript. 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. Data Availability Statement: The data presented in this study are contained within the article. Conflicts of Interest: The authors declare no conflict of interest. References 1. Furman, D.; Campisi, J.; Verdin, E.; Carrera-Bastos, P.; Targ, S.; Franceschi, C.; Ferrucci, L.; Gilroy, D.W.; Fasano, A.; Miller, G.W.; et al. Chronic Inflammation in the Etiology of Disease across the Life Span. Nat. Med. 2019, 25, 1822–1832. [CrossRef] [PubMed] 2. Serhan, C.N.; Chiang, N.; van Dyke, T.E. Resolving Inflammation: Dual Anti-Inflammatory and pro-Resolution Lipid Mediators. Nat. Rev. Immunol. 2008, 8, 349–361. [CrossRef] [PubMed] 3. Ayoub, S.S. Fundamentals of Inflammation; Serhan, C.N., Ward, P.A., Gilroy, D.W., Eds.; Cambridge University Press: Cambridge, UK, 2010; ISBN 9781139195737. 4. Buckley, C.D.; Gilroy, D.W.; Serhan, C.N. Proresolving Lipid Mediators and Mechanisms in the Resolution of Acute Inflammation. Immunity 2014, 40, 315–327. [CrossRef] 5. Serhan, C.N. Discovery of Specialized Pro-Resolving Mediators Marks the Dawn of Resolution Physiology and Pharmacology. Mol. Asp. Med. 2017, 58, 1–11. [CrossRef] [PubMed] 6. Dalli, J. Does Promoting Resolution Instead of Inhibiting Inflammation Represent the New Paradigm in Treating Infections? Mol. Asp. Med. 2017, 58, 12–20. [CrossRef] 7. Serhan, C.N. Treating Inflammation and Infection in the 21st Century: New Hints from Decoding Resolution Mediators and Mechanisms. FASEB J. 2017, 31, 1273–1288. [CrossRef] [PubMed] 8. Wan, M.; Tang, X.; Stsiapanava, A.; Haeggström, J.Z. Biosynthesis of Leukotriene B4. Semin. Immunol. 2017, 33, 3–15. [CrossRef] [PubMed] 9. Mancini, J.A.; Abramovitz, M.; Cox, M.E.; Wong, E.; Charleson, S.; Perrier, H.; Wang, Z.; Prasit, P.; Vickers, P.J. 5-Lipoxygenase- Activating Protein Is an Arachidonate Binding Protein. FEBS Lett. 1993, 318, 277–281. [CrossRef] 10. Peters-Golden, M.; Brock, T.G. 5-Lipoxygenase and FLAP. Prostaglandins Leukot. Essent. Fat. Acids 2003, 69, 99–109. [CrossRef] 11. Evans, J.F.; Ferguson, A.D.; Mosley, R.T.; Hutchinson, J.H. What’s All the FLAP about?: 5-Lipoxygenase-Activating Protein Inhibitors for Inflammatory Diseases. Trends Pharmacol. Sci. 2008, 29, 72–78. [CrossRef] 12. Gerstmeier, J.; Weinigel, C.; Rummler, S.; Rådmark, O.; Werz, O.; Garscha, U. Time-resolved in Situ Assembly of the Leukotriene- synthetic 5-lipoxygenase/5-lipoxygenase-activating Protein Complex in Blood Leukocytes. FASEB J. 2016, 30, 276–285. [CrossRef] 13. Abramovitz, M.; Wong, E.; Cox, M.E.; Richardson, C.D.; Li, C.; Vickers, P.J. 5-Lipoxygenase-Activating Protein Stimulates the Utilization of Arachidonic Acid by 5-Lipoxygenase. Eur. J. Biochem. 1993, 215, 105–111. [CrossRef] Int. J. Mol. Sci. 2022, 23, 3140 22 of 23 14. Basavarajappa, D.; Wan, M.; Lukic, A.; Steinhilber, D.; Samuelsson, B.; Radmark, O. Roles of Coactosin-like Protein (CLP) and 5-Lipoxygenase-Activating Protein (FLAP) in Cellular Leukotriene Biosynthesis. Proc. Natl. Acad. Sci. USA 2014, 111, 11371–11376. [CrossRef] 15. Rudberg, P.C.; Tholander, F.; Thunnissen, M.M.G.M.; Samuelsson, B.; Haeggstrom, J.Z. Leukotriene A4 Hydrolase: Selective Abrogation of Leukotriene B4 Formation by Mutation of Aspartic Acid 375. Proc. Natl. Acad. Sci. USA 2002, 99, 4215–4220. [CrossRef] 16. Thunnissen, M.M.G.M.; Nordlund, P.; Haeggström, J.Z. Crystal Structure of Human Leukotriene A4 Hydrolase, a Bifunctional Enzyme in Inflammation. Nat. Struct. Biol. 2001, 8, 131–135. [CrossRef] 17. Rinaldo-Matthis, A.; Haeggström, J.Z. Leukotriene A4 Hydrolase and Leukotriene C4 Synthase. In Lipoxygenases in Inflammation; Springer: Frankfurt, Germany, 2016; pp. 31–46. [CrossRef] 18. Stsiapanava, A.; Samuelsson, B.; Haeggström, J.Z. Capturing LTA4 Hydrolase in Action: Insights to the Chemistry and Dynamics of Chemotactic LTB4 Synthesis. Proc. Natl. Acad. Sci. USA 2017, 114, 9689–9694. [CrossRef] 19. Rudberg, P.C.; Tholander, F.; Thunnissen, M.M.G.M.; Haeggström, J.Z. Leukotriene A4 Hydrolase/Aminopeptidase: Glutamate 271 is a catalytic residue with specific roles in two distinct enzyme mechanisms. J. Biol. Chem. 2002, 277, 1398–1404. [CrossRef] 20. Rudberg, P.C.; Tholander, F.; Andberg, M.; Thunnissen, M.M.G.M.; Haeggström, J.Z. Leukotriene A4 Hydrolase: Identification of a common carboxylate recognition site for the epoxide hydrolase and aminopeptidase substrates. J. Biol. Chem. 2004, 279, 27376–27382. [CrossRef] 21. Mueller, M.J.; Blomster, M.; Oppermann, U.C.; Jornvall, H.; Samuelsson, B.; Haeggstrom, J.Z. Leukotriene A4 Hydrolase: Protection from Mechanism-Based Inactivation by Mutation of Tyrosine-378. Proc. Natl. Acad. Sci. USA 1996, 93, 5931–5935. [CrossRef] 22. Mueller, M.J.; Blomster Andberg, M.; Samuelsson, B.; Haeggstrom, J.Z. Leukotriene A4 Hydrolase Mutation of Tyrosine 378 Allows Conversion of Leukotriene A4 into an Isomer of Leukotriene B4. J. Biol. Chem. 1996, 271, 24345–24348. [CrossRef] 23. Mu, X.; Xu, D. QM/MM Molecular Dynamics Investigations of the Substrate Binding of Leucotriene A4 Hydrolase: Implication for the Catalytic Mechanism. J. Phys. Chem. B 2018, 122, 7253–7263. [CrossRef] 24. Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260–7268. [CrossRef] 25. Elstner, M.; Frauenheim, T.; Kaxiras, E.; Seifert, G.; Suhai, S. A Self-Consistent Charge Density-Functional Based Tight-Binding Scheme for Large Biomolecules. Phys. Status Solidi 2000, 217, 357–376. [CrossRef] 26. Elstner, M. The SCC-DFTB Method and Its Application to Biological Systems. Theor. Chim. Acta 2006, 116, 316–325. [CrossRef] 27. Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M. A QM/MM Implementation of the Self-Consistent Charge Density Functional Tight Binding (SCC-DFTB) Method. J. Phys. Chem. B 2001, 105, 569–585. [CrossRef] 28. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera? A Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605–1612. [CrossRef] [PubMed] 29. Dijkstra, E.W. Invariance and Non-Determinacy. Numer. Math. 1959, 1, 269–271. [CrossRef] 30. Olsson, M.H.M.; Søndergaard, C.R.; Rostkowski, M.; Jensen, J.H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput. 2011, 7, 525–537. [CrossRef] 31. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. [CrossRef] 32. Sherwood, P.; de Vries, A.H.; Guest, M.F.; Schreckenbach, G.; Catlow, C.R.A.; French, S.A.; Sokol, A.A.; Bromley, S.T.; Thiel, W.; Turner, A.J.; et al. QUASI: A General Purpose Implementation of the QM/MM Approach and Its Application to Problems in Catalysis. J. Mol. Struct. THEOCHEM 2003, 632, 1–28. [CrossRef] 33. Metz, S.; Kästner, J.; Sokol, A.A.; Keal, T.W.; Sherwood, P. ChemShell—A Modular Software Package for QM/MM Simulations. WIREs Comput. Mol. Sci. 2014, 4, 101–110. [CrossRef] 34. Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165–169. [CrossRef] 35. Case, D.A.; Ben-Shalom, I.Y.; Brozell, S.R.; Cerutti, D.S.; Cheatham, T.E., III; Cruzeiro, V.W.D.; Darden, T.A.; Duke, R.E.; Ghoreishi, D.; Gilson, M.K.; et al. Amber 18; University of California: San Francisco, CA, USA, 2018. 36. Smith, W.; Forester, T.R. DL_POLY_2.0: A General-Purpose Parallel Molecular Dynamics Simulation Package. J. Mol. Graph. 1996, 14, 136–141. [CrossRef] 37. Maier, J.A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K.E.; Simmerling, C. Ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from Ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. [CrossRef] 38. Wang, J.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157–1174. [CrossRef] 39. Bayly, C.I.; Cieplak, P.; Cornell, W.D.; Kollman, P.A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269–10280. [CrossRef] 40. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.; Nakatsuji, H.; et al. Gaussian09; Gaussian, Inc.: Wallingford, CT, USA, 2009. Int. J. Mol. Sci. 2022, 23, 3140 23 of 23 41. Bakowies, D.; Thiel, W. Hybrid Models for Combined Quantum Mechanical and Molecular Mechanical Approaches. J. Phys. Chem. 1996, 100, 10580–10594. [CrossRef] 42. Senn, H.M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int. Ed. 2009, 48, 1198–1229. [CrossRef] 43. De Vries, A.H.; Sherwood, P.; Collins, S.J.; Rigby, A.M.; Rigutto, M.; Kramer, G.J. Zeolite Structure and Reactivity by Combined Quantum-Chemical−Classical Calculations. J. Phys. Chem. B 1999, 103, 6133–6141. [CrossRef] 44. Kästner, J.; Thiel, S.; Senn, H.M.; Sherwood, P.; Thiel, W. Exploiting QM/MM Capabilities in Geometry Optimization: A Microiterative Approach Using Electrostatic Embedding. J. Chem. Theory Comput. 2007, 3, 1064–1072. [CrossRef] 45. Billeter, S.R.; Turner, A.J.; Thiel, W. Linear Scaling Geometry Optimisation and Transition State Search in Hybrid Delocalised Internal Coordinates. Phys. Chem. Chem. Phys. 2000, 2, 2177–2186. [CrossRef] 46. Liu, D.C.; Nocedal, J. On the Limited Memory BFGS Method for Large Scale Optimization. Math. Program. 1989, 45, 503–528. [CrossRef] 47. Nocedal, J. Updating Quasi-Newton Matrices with Limited Storage. Math. Comput. 1980, 35, 773–782. [CrossRef] 48. Baker, J. An Algorithm for the Location of Transition States. J. Comput. Chem. 1986, 7, 385–395. [CrossRef] 49. Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. Search for Stationary Points on Surfaces. J. Phys. Chem. 1985, 89, 52–57. [CrossRef] 50. Kästner, J.; Carr, J.M.; Keal, T.W.; Thiel, W.; Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. J. Phys. Chem. A 2009, 113, 11856–11865. [CrossRef] [PubMed] 51. Becke, A.D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. [CrossRef] 52. Hariharan, P.C.; Pople, J.A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213–222. [CrossRef] 53. Schautz, F.; Flad, H.-J.; Dolg, M. Quantum Monte Carlo Study of Be 2 and Group 12 Dimers M 2 (M = Zn, Cd, Hg). Theor. Chim. Acta 1998, 99, 231–240. [CrossRef] 54. Case, D.A.; Ben-Shalom, I.Y.; Brozell, S.R.; Cerutti, D.S.; Cheatham, T.E., III; Cruzeiro, V.W.D.; Darden, T.A.; Duke, R.E.; Ghoreishi, D.; Giambasu, G.; et al. Amber19; University of California: San Francisco, CA, USA, 2019. 55. Seminario, J.M. Calculation of Intramolecular Force Fields from Second-Derivative Tensors. Int. J. Quantum Chem. 1996, 60, 1271–1277. [CrossRef] 56. Li, P.; Merz, K.M. MCPB.Py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016, 56, 599–604. [CrossRef] [PubMed] 57. Le Grand, S.; Götz, A.W.; Walker, R.C. SPFP: Speed without Compromise—A Mixed Precision Model for GPU Accelerated Molecular Dynamics Simulations. Comput. Phys. Commun. 2013, 184, 374–380. [CrossRef] 58. Götz, A.W.; Williamson, M.J.; Xu, D.; Poole, D.; le Grand, S.; Walker, R.C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8, 1542–1555. [CrossRef] [PubMed] 59. Ryckaert, J.P.; Ciccotti, G.; Berendsen, H.J.C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327–341. [CrossRef] 60. Leach, A.R. Molecular Modelling: Principles and Applications; Pearson Education: London, UK, 2001. 61. Berendsen, H.J.C.; Postma, J.P.M.; van Gunsteren, W.F.; Dinola, A.; Haak, J.R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. [CrossRef]