International Journal of Molecular Sciences Article Hydroperoxidation of Docosahexaenoic Acid by Human ALOX12 and pigALOX15-mini-LOX 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, Bellaterra, 08193 Barcelona, Spain 2 Arquebio SL, Carrer de Álava 51, 08005 Barcelona, Spain 3 Institut de Biotecnologia i Biomedicina (IBB), Universitat Autònoma de Barcelona, Bellaterra, 08193 Barcelona, Spain * Correspondence: angels.gonzalez@uab.cat Abstract: Human lipoxygenase 12 (hALOX12) catalyzes the conversion of docosahexaenoic acid (DHA) into mainly 14S-hydroperoxy-4Z,7Z,10Z,12E,16Z,19Z-docosahexaenoic acid (14S-H(p)DHA). This hydroperoxidation reaction is followed by an epoxidation and hydrolysis process that finally leads to maresin 1 (MaR1), a potent bioactive specialized pro-resolving mediator (SPM) in chronic inflammation resolution. By combining docking, molecular dynamics simulations, and quantum mechanics/molecular mechanics calculations, we have computed the potential energy profile of DHA hydroperoxidation in the active site of hALOX12. Our results describe the structural evolution of the molecular system at each step of this catalytic reaction pathway. Noteworthy, the required stereospecificity of the reaction leading to MaR1 is explained by the configurations adopted by DHA bound to hALOX12, along with the stereochemistry of the pentadienyl radical formed after the first step of the mechanism. In pig lipoxygenase 15 (pigALOX15-mini-LOX), our calculations suggest that 14S-H(p)DHA can be formed, but with a stereochemistry that is inadequate for MaR1 biosynthesis. Keywords: hydroperoxidation mechanism; human platelet ALOX12; enzyme catalysis; molecular dynamics simulations; QM/MM calculations Citation: Canyelles-Niño, M.; 1. Introduction González-Lafont, À.; Lluch, J.M. Mammalian lipoxygenases (LOXs) form a heterogeneous family of non-heme iron- Hydroperoxidation of containing enzymes that catalyze the dioxygenation of polyunsaturated fatty acids (PUFAs) Docosahexaenoic Acid by Human with cis-methylene-interrupted double bonds. The peroxidation process is highly regio- ALOX12 and pigALOX15-mini-LOX. and stereospecific, leading to the corresponding conjugated hydroperoxides [1]. Then, Int. J. Mol. Sci. 2023, 24, 6064. along the biological pathway, those hydroperoxy fatty acids formed by the different LOXs https://doi.org/10.3390/ijms24076064 isoforms generate a wide variety of bioactive lipid mediators that are involved in a plethora Academic Editor: Jan Korabecny of pro- and anti-inflammatory processes [1–3]. The first lipoxygenase discovered in animal tissues by Nobel Laureate Samuelsson Received: 13 February 2023 Revised: 15 March 2023 and his colleague Hamberg was the human 12-lipoxygenase of blood platelets [4]. This Accepted: 17 March 2023 enzyme was named 12S-LOX because it oxidizes arachidonic acid (AA) to generate almost Published: 23 March 2023 exclusively one regio- and stereospecific isomer, the 12S-hydro(per)oxyeicosa-5,8,10,14- tetraenoic acid (12S-H(p)ETE) [5]. However, this nomenclature was confusing because the regio- and stereospecificity of LOX isoforms depend on the substrate and are species- specific [6,7]. Moreover, this nomenclature is not related to the biological function of LOX Copyright: © 2023 by the authors. isoforms. For this reason, terminology based on the gene encoding each LOX isoform is Licensee MDPI, Basel, Switzerland. more recommended. In humans, the platelet-type 12-LOX (12S-LOX) is encoded by the This article is an open access article ALOX12 gene and is named human ALOX12 (hALOX12) [7]. The hALOX12 peroxidation distributed under the terms and reaction is initiated by H-abstraction of the H10 atom of AA by the Fe(III)-OH− cofactor. conditions of the Creative Commons Next, an oxygen molecule inserts on C12 following an antarafacial approach and gives a Attribution (CC BY) license (https:// hydroperoxide with S stereochemistry. Finally, the reduced product, 12S-HETE, is obtained creativecommons.org/licenses/by/ by the action of glutathione peroxidase. 4.0/). Int. J. Mol. Sci. 2023, 24, 6064. https://doi.org/10.3390/ijms24076064 https://www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2023, 24, 6064 2 of 29 Human ALOX12 is the unique LOX isoform expressed in human platelets [4]. Nowa- days this LOX has been found in platelets of several other mammals such as bovine, mouse, rat, rabbit, sheep, and dog [8–11] but is not expressed in porcine platelets [12]. Through hALOX12 oxidized metabolites, this enzyme regulates platelet activity in hemostasis and thrombosis, and controls platelet role in inflammatory processes. According to several studies with human patients, hALOX12 and its derived oxylipin 12S-HETE seem to present prothrombotic effects [13]. However, the prothrombotic role of 12S-HETE in vivo is not fully understood yet. There are contradictory experimental results with animal models in which 12S-HETE could either promote or inhibit or even not participate in platelet aggregation [14,15]. In this respect, the first selective hALOX12 inhibitor, ML355, has been shown to impair thrombus formation and vessel occlusion in vivo with minimal effects on hemostasis [16,17]. It has also been recognized that platelets are involved in the pathogenesis of several diseases, such as asthma, cancer, neurodegenerative disorders, gastrointestinal and hepatic inflammation, insulin resistance, and atherosclerosis, in which chronic inflammation plays an underlying role [7,18]. So, 12-HETE has been described as a pro-inflammatory molecule. However, it is not always clear whether the 12-HETE present in the inflamed tissue originated from platelets or from other cells because hALOX12 is also expressed in leucocytes, human islets, and epidermal keratinocytes [19]. On the other hand, hALOX12, in conjunction with other LOX isoforms, is involved in the transcellular biosynthesis of AA-derived anti-inflammatory agents such as lipoxins [20]. Human ALOX12 also catalyzes the oxygenation of the ω-3 fatty acid 4,7,10,13,16,19- docosahexaenoic acid (DHA). DHA can be converted by hALOX12 into two different products, depending on the experimental conditions. In this respect, there are different results in the literature concerning the regioselectivity of the DHA oxygenation reaction with hALOX12. According to Kühn and co-workers [21], the ω-9 hydroperoxide prod- uct 14-hydroperoxydocosahexaenoic acid (14-H(p)DHA) is dominant over DHA. Those experimental results did not show significant amounts of any other product with this enzyme. Thus, hALOX12 catalysis with DHA was classified as singular ω-9 oxygenase activity. Two years later, Holman and co-workers [22] published new kinetic results that confirmed 14-H(p)DHA as the main product of DHA peroxidation by hALOX12. However, the obtained percentage of 11-HDHA:14-HDHA (29:71) products (the reduced forms of 11-hydroperoxydocosahexaenoic acid (11-H(p)DHA) and 14-H(p)DHA, respectively) was higher than previously reported. According to the authors, those differences appear be- cause they expressed hALOX12 in insect cells, whereas Kühn and co-workers used E. coli as the expression system [21]. In a following paper, Holman and co-workers report a smaller in vitro proportion of 11S-HDHA (19%) to 14S-HDHA (81%), and an even smaller ex vivo percentage (6:94) [23]. The DHA peroxidation mechanism is proposed to initiate with an H-abstraction reaction (described as a proton-coupled electron transfer, PCET [24,25]) of H12 or H9, in this case, by the Fe(III)-OH− cofactor. Next, an oxygen molecule inserts on C14 or C11 of DHA, following an antarafacial approach, and gives the corresponding hydroperoxides with S stereochemistry, which are immediately reduced to 14-HDHA and 11-HDHA, respectively. Moreover, steady-state kinetic measurements of this catalytic mechanism of hALOX12 have demonstrated that the rates of substrate capture and product release are similar between DHA and AA [23,26]. The biological relevance of DHA-derived oxylipins, 11-HDHA and 14-HDHA, is cur- rently a focus of interest. The cardiovascular benefits of a diet rich in ω-3 fatty acids such as DHA have been known for decades, but the mechanisms involved in their anti-aggregatory effects on platelets are still under study. Recent in vitro and ex vivo experiments [27] with mice and human-washed platelets have confirmed that DHA, as well as 11-HDHA and 14-HDHA, inhibit collagen-induced platelet aggregation. Interestingly, the in vivo thrombosis data in mice showed that only acute concentrations of 11-HDHA or 14-HDHA attenuated overall thrombus formation. This fact suggests that under physiological con- ditions (at similar basal ALOX12-derived oxylipin levels than those used in the assays mentioned), 11-HDHA and 14-HDHA may play a central role in the regulation of throm- Int. J. Mol. Sci. 2023, 24, 6064 3 of 29 bus formation by activating Protein Kinase A (PKA). Nevertheless, the authors of the experimental study do not discard additional complex mechanisms driven by special- ized pro-resolving mediators (SPMs), also derived from DHA, that could contribute to antithrombotic effects. Maresin 1 (7R,14S-dihydroxy-4Z,8E,10E,12Z,16Z,19Z- docosahex- aenoic acid, MaR1) is an SPM derived from 14S-H(p)DHA that is subsequently converted to 13S,14S-epoxy-4Z,7Z,9E,11E,16Z,19Z-docosahexaenoic acid through an epoxidation mech- anism. Then, this epoxide intermediate is proposed to be enzymatically hydrolyzed to form the stereochemistry of bioactive MaR1 [23,28]. MaR1 is a potent SPM that actively participates in the inflammation-resolution phase of many different diseases, and it has also been described as a novel antiplatelet agent crucial in the resolution of inflammation in cardiovascular injuries [29,30]. At the molecular level, the structural basis of hALOX12’s product selectivity has been much less studied than in the case of rabbit and human ALOX15. Remarkably, some muta- genesis experiments back in 2010 already showed that the triad concept (which established that the three residues, Phe353, Ile418, and Phe593, at the bottom of rabbit ALOX15’s cavity were critical for AA positioning) only partially applies to hALOX12 [31,32]. Instead, this hypothesis has been supported several times for ALOX15 orthologs. By mutagenesis of the residues of the triad, an AA 15-lipoxigenating ALOX15 ortholog has been converted into an AA 12-lipoxygenating enzyme or vice versa. However, a U-shaped binding pocket [33–35] was introduced to complete the description of the substrate binding cavity of LOX isoforms versus the triad concept [6]. Moreover, the regioselectivity of LOX isoforms has also been shown to depend on the orientation of the substrate bound to the enzyme. That is, which end of the substrate (carboxylate-head or hydrophobic-tail) is innermost in the cavity [6,35]. Recently, Holman and co-workers [22] have thoroughly studied, by site-directed mutagene- sis and, using docking calculations, the main interactions of AA with hALOX12 active-site residues, establishing the main differences between hALOX12 and hALOX15 regioselectivities. In that work, the double mutant A417I/V418M of hALOX12 increased 15-HpETE production by only 24 ± 2%, whereas the single mutation of hALOX15 (I418A) augmented the yield of 12-HETE to 94% [31]. Also, the docking simulations by Holman and co-workers supported the U-shaped binding mode [35] for AA in hALOX12 with a tail-first orientation [22]. How- ever, the calculations in that paper were performed using a model of the catalytic domain of hALOX12 built from the crystallographic coordinates of an N-terminal truncation variant of pigALOX15 [36], so-called pigALOX15-mini-LOX. As the same authors explain, the X-ray structure for hALOX12‘s catalytic domain, available at the Protein Data Bank (PDB code 3D3L), is not adequate for modeling the AA/hALOX12 complex. The porcine LOX isoform is encoded by the ALOX15 gene, and for that reason, it is called pigALOX15. It presents a high degree (86%) of sequence conservation with ra- bALOX15, but much less similarity with hALOX12 (66%). This porcine LOX was purified from the cytosolic fraction of leukocytes [37], and it was later also mainly found in the anterior pituitary gland of pigs [38]. As with the other ALOX15 orthologs, pigALOX15 presents dual specificity when reacting with AA. The main pigALOX15 peroxidation prod- ucts are 12S- and 15S-HpETE in a ratio of about 10:1 [39,40]. However, as mentioned above, this regioselectivity has been inverted by site-directed mutagenesis based on the triad concept. In a study by some of the authors [41], the complete wild-type pigALOX15 and its pigALOX15-mini-LOX variant were expressed, and both showed the same regiospecificity, mainly producing 12-H(p)ETE. However, catalytic activity was strongly reduced for the isolated catalytic domain. The specificity switch of pigALOX15-mini-LOX induced by experimental mutagenesis of Val418 and Val419 sequence determinants was demonstrated. Moreover, docking, MD simulations, and quantum mechanics/molecular mechanics calcu- lations using a solvated model of wild-type and the Val418Ile + Val419Met double mutant of pigALOX15-mini-LOX were performed. In agreement with experiments, the lower energy barriers in the wild-type LOX were found for the H10 -abstraction process leading to 12-lipoxygenation whereas the double mutant was dominantly 15-lipoxygenating with lower barriers for the H13 -abstraction. Int. J. Mol. Sci. 2023, 24, 6064 4 of 29 In this study, we have built for the first time an in silico model of the DHA/hALOX12 complex using the AlphaFold [42,43] server to retrieve the initial coordinates of the human LOX protein. Docking and MD simulations have been performed to simulate the binding mode of DHA at the active site of hALOX12. Next, we have calculated within a QM/MM scheme the potential energy profile of the hydroperoxidation mechanistic steps that convert DHA to the fully characterized 14S-hydroperoxy-4Z,7Z,10Z,12E,16Z,19Z-docosahexaenoic acid (14S-H(p)DHA). Special attention has been devoted to the stereochemistry of the reaction because the biosynthesis of maresin 1, which is the potent bioactive SPM derived from 14S-H(p)DHA, is highly stereospecific. For the sake of comparison, the viability of DHA peroxidation by pigALOX15 has also been analyzed using the crystallographic coordinates of an N-terminal truncation variant of pigALOX15 (pigALOX15-mini-LOX) that has been previously used as a model of the hALOX12 structure. 2. Results and Discussion Docking calculations were carried out for the DHA/hALOX12 and DHA/pigALOX15- mini-LOX complexes. The binding mode of the substrate with the best score in the active site of each LOX was selected to initiate the corresponding MD simulations. 2.1. MD Simulations of the DHA/hALOX12 Complex In this section, the achieved outcomes from the two MD trajectories of the solvated DHA/hALOX12 complex are presented. In Figure 1, we have plotted the protein back- bone RMSDs versus time. The protein RMSD along the first replica remains quite stable during 170 ns. Then, a conformational change occurs that increases the RMSD. This con- formational change corresponds to an opening of the PLAT (Polycystin-1,Lipoxygenase, Alpha-Toxin) domain (N-terminal domain from residues 1–110) away from the catalytic domain (C-terminal domain) (Figure 2). The average distance between the geometrical centers of the two domains changes from 42.3 Å to 47.2 Å when the domains separate (see Figure S1). The angle between the longest axes of the two ellipsoids that enclose the PLAT and the catalytic domain, respectively, is 38.2◦ in the closed conformation and 58.8◦ in the open structure. Experimentally, the flexibility of the PLAT domain of hALOX12, described as a pendulum-like movement, was proposed for fitting SAXS measurements [44,45]. The C-skeleton of DHA also experiences a sudden conformational change just before the N- terminal displacement (see Figure 1, replica 1). In the second MD replica, the PLAT domain remains more tightly bound to the catalytic domain in a close conformation (see Figure 1, replica 2). Moreover, in this second MD replica, the DHA substrate shows the same binding mode all along the trajectory. In the hALOX12 active site, DHA adopts a U-shaped binding mode [35], like that described for AA by Holman and co-workers using docking calculations [22], while the two-domain LOX remains closed. There is a modification of the DHA pose bound to hALOX12 when the PLAT domain separates. An overlay of both DHA binding modes is depicted in Figure 3. In Figure 4, the interatomic distances of the main interactions between the carboxylate head of DHA and hALOX12 residues have been plotted along the two MD trajectories. As DHA has two more carbons than AA, its carboxylate group is hydrogen-bonded to Arg402 instead of being bound to His596 like AA (see Figure 3). Arg402 is a more solvent-exposed residue than His596 that can only be reached by DHA being a longer substrate. The average hydrogen-bond distance between Arg402 and the carboxylate of DHA (2.41 Å and 2.12 Å in replicas 1 and 2, respectively) is similar in both MD replicas, but in replica 2, the interaction remains more stable all along the trajectory. The carboxylate group of DHA also interacts with Gln406, but this interaction is less favorable and more fluctuating than with Arg402. In turn, Arg402 and Gln406 interact by hydrogen bonds along the MD trajectory, as can be seen in Figure 5 for replica 1. However, Gln406 flips repeatedly away from the cavity, losing contact with DHA’s carboxylate and Arg402. In contrast, Arg402 remains more fixed along the trajectory due to its electro- Int. J. Mol. Sci. 2023, 24, 6064 5 of 29 Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 5 of 32 static interactions with Glu175 of the α2-helix. Several water molecules also contribute to stabilizing the substrate’s carboxylate head. Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 5 of 32 Protein FigureFigure 1. Protein andsubstrate and 1. Protein substrate and backbone backbone substrate RMSDs RMSDs backbone versus RMSDsversus timefor versustime time for for the the the MD MDMD replica replica replica 1 replica and 1 and 1 and replica replica 2 2 in2 in hALOX12. in hALOX12. hALOX12. 2. PLAT FigureFigure domain 2. PLAT of hALOX12 domain of hALOX12in theinclosed (in green) the closed and and (in green) in the in open (in orange) the open conformation. (in orange) confor- mation. The The C-terminal C-terminal domain domain is shown in isgrey. shown in grey. Figure 2. PLAT domain of hALOX12 in the closed (in green) and in the open (in orange) confor- mation. The C-terminal domain is shown in grey. In the hALOX12 active site, DHA adopts a U-shaped binding mode [35], like that described for AA by Holman and co-workers using docking calculations [22], while the two-domain LOX remains closed. There is a modification of the DHA pose bound to hA- Int. J. Mol. Sci. 2023, 24, 6064 6 of 29 LOX12 when the PLAT domain separates. An overlay of both DHA binding modes is de- picted in Figure 3. Figure 3. DHA binding mode in the closed (in green) and the open (in red) conformations of hA- Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW Figure 3. DHA binding mode in the closed (in green) and the open (in red) conformations of7 of 32 LOX12. The interaction between Arg402 and the carboxylate group of DHA is shown as well as the hALOX12. The interaction between Arg402 and the carboxylate group of DHA is shown as well as rotation of Phe414. the rotation of Phe414. In Figure 4, the interatomic distances of the main interactions between the carbox- ylate head of DHA and hALOX12 residues have been plotted along the two MD trajecto- ries. As DHA has two more carbons than AA, its carboxylate group is hydrogen-bonded to Arg402 instead of being bound to His596 like AA (see Figure 3). Arg402 is a more sol- vent-exposed residue than His596 that can only be reached by DHA being a longer sub- strate. The average hydrogen-bond distance between Arg402 and the carboxylate of DHA (2.41 Å and 2.12 Å in replicas 1 and 2, respectively) is similar in both MD replicas, but in replica 2, the interaction remains more stable all along the trajectory. Figure Figure 4. Distancesofofa acarboxylate’s 4. Distances carboxylate’s oxygen oxygenatom atomofofDHA DHA to to thethe closest hydrogen closest atomatom hydrogen in Arg402 in Arg402 and Gln406 versus time for the MD replica 1 and replica 2 in hALOX12. and Gln406 versus time for the MD replica 1 and replica 2 in hALOX12. The carboxylate group of DHA also interacts with Gln406, but this interaction is less favorable and more fluctuating than with Arg402. In turn, Arg402 and Gln406 interact by hydrogen bonds along the MD trajectory, as can be seen in Figure 5 for replica 1. However, Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 8 of 32 Int. J. Mol. Sci. 2023, 24, 6064 7 of 29 Figure Figure 5. 5. Distancebetween Distance betweenthe the oxygen oxygen atom atomininGln406 Gln406and thethe and closest hydrogen closest atomatom hydrogen in Arg402 in Arg402 versus versus time time forfor theMD the MDreplica replica 11 in in hALOX12. hALOX12. Holman and co-workers mutated the Arg402 residue by a hydrophobic Leu and did notHolman observe and co-workers any significant mutated change in thethe Arg402ofresidue reactivity hALOX12 bywith a hydrophobic Leu and did DHA [22]. However, notthis observe resultanydoessignificant not necessarilychange meanin the thatreactivity Arg402 does of hALOX12 not interactwithwithDHA [22]. However, the carboxylate thisgroup resultofdoes not necessarily mean that Arg402 does not DHA in WT hALOX12, as was suggested by the authors of the mutagenesis interact with the carboxylate group of DHA In experiment. in the WTArg402Leu hALOX12,mutant, as wasthe suggested by the authors DHA carboxylate group of the reinforce might mutagenesis its ex- interaction periment. with In the Gln406 andmutant, Arg402Leu the closerthewater DHAmolecules. carboxylate Thus, groupthe might changereinforce of hydrogen- its inter- bonded action withpartner Gln406(mainly and theArg402 closerin WT hALOX12 water molecules.byThus, Gln406theinchange the Arg402Leu mutant) of hydrogen-bonded would not imply any relocation of DHA. partner (mainly Arg402 in WT hALOX12 by Gln406 in the Arg402Leu mutant) would not DHA also establishes π–π interactions between its six double bonds and aromatic imply any relocation of DHA. residues of the hALOX12 active site that stabilize the substrate binding mode. In Figure 6, weDHAplotted also theestablishes distances between π–π interactions Phe174 and the between its six atom closer carbon double of ∆bonds 4 and ∆and aromatic 7 double residues bonds of the the along hALOX12 active site MD trajectory that stabilize of replica 1. This π–πthe substrate binding the stacking stabilizes mode. In Figure 6, ∆4 double webond plotted the distances between Phe174 and the closer carbon atom during the first 150 ns of the trajectory until the DHA substrate moves to the cavity of Δ 4 and Δ7 double bonds alongThen, entrance. the MDPhe174trajectory switches ofits replica stacking 1. This π–π stacking interactions ∆4 to the ∆the from thestabilizes 7 double Δ4 double bond.bond In Figure during S2, the the first 150restnsofofthethestacking trajectoryinteractions until thebetween the substrate DHA substrate movesand tothethe enzyme cavity en- residues are plotted. As can be observed, there are also π–π interactions trance. Then, Phe174 switches its stacking interactions from the Δ to the Δ double bond. 4 between 7 Phe352 In Figure ∆19 the and the S2, double restbondof theandstacking His365 and the ∆10 double interactions bond.the between Thesubstrate stacking interaction and the enzyme of DHA with Phe414, which has been described for AA, is established with two double bonds residues are plotted. As can be observed, there are also π–π interactions between Phe352 (∆13 and19 ∆16 ). In the MD replica 1, all those interactions weaken when the PLAT domain andseparates the Δ double bond and His365 and the Δ10 double bond. The stacking interaction of from the catalytic domain. The principal hALOX12 residues interacting with DHA DHA’swith Phe414, double bonds which has been by stacking described interactions are for AA,inisFigure plotted established 7 for thewith two double two binding bonds modes (Δobserved 13 and Δ along 16 ). In the theMD MDreplica replica 1, replica 1. In all those interactions 2, the weaken with stacking interactions when the PLAT Phe174 domain remain separates stable afterfrom80 nsthe(Figure catalytic domain. 6). In Figure The S2, itprincipal hALOX12 can be observed that residues interacting the interactions with with Phe352double DHA’s and Phe414bonds remain very stable by stacking all along the interactions aretrajectory. plotted in In Figure contrast,7 the forπ–πthe two stacking with His365 is not present in this MD replica. binding modes observed along the MD replica 1. In replica 2, the stacking interactions The PLAT with Phe174 movement remain stable also aftercorrelates 80 ns (Figurewith structural 6). In Figurechanges S2, in somebeα-helices it can observed of the that the catalytic domain. In this respect, it is noteworthy that there is a loss of the secondary interactions with Phe352 and Phe414 remain very stable all along the trajectory. In con- structure of the α-helix from residue 408 to residue 418, which flips the sidechain of Phe414. trast, the π–π stacking with His365 is not present in this MD replica. With the displacement of the PLAT domain, Phe414 rotates, blocking the U-shaped cavity bottomPLAT The movement and opening a new also cavitycorrelates to accommodate with structural the DHA tail. changes in some This rotation of α-helices the Phe414of the catalytic sidechaindomain. In thisinrespect, is depicted Figures it3,is7 noteworthy and S3. Then, that there Ile408 is a loss ofthe approaches theDHA secondary tail thatstruc- ture of the α-helix simultaneously fromaround rotates residue the408 to 18residue C17 –C 418, which single bond, losing theflips the sidechain U-shaped of Phe414. binding mode. With the displacement of the PLAT domain, Phe414 rotates, blocking the U-shaped cavity bottom and opening a new cavity to accommodate the DHA tail. This rotation of the Phe414 sidechain is depicted in Figures 3, 7, and S3. Then, Ile408 approaches the DHA tail that simultaneously rotates around the C17–C18 single bond, losing the U-shaped binding mode. Int.Int. J. Mol. Sci. J. Mol. 2023,24, Sci.2023, 6064 24,Int. J. Mol. x FOR PEERSci.REVIEW 2023, 24, x FOR PEER REVIEW 9 of 32 8 of 29 9 of 32 Distances Figure6.6.Distances Figure between Figure between Phe174 6.Phe174 Distances side chainchain sidebetween Phe174 and and the the chain side closest closestofatom atomand of theΔclosest the 4 and 4∆ theΔatomand ∆bonds 7 of the 7 double double Δ4 and Δbonds 7 double bonds versustime versus timeforfor versus thethe MD MD time1 for replica replica 1the and andMD replica replica replica 2 1inand replica 2 in hALOX12. hALOX12. 2 in hALOX12. Figure7.7.Main Figure Figure Mainstacking stacking 7. Main stacking interactions interactions between interactions between hALOX12 hALOX12 between hALOX12 residues residues and residues and DHA. DHA. Theand DHA.residues protein The protein residuesThe protein and residues and DHA are shown and in DHA green are for shown the in closed green for the conformation, closed and conformation, in orange and and red, in orange and respectively, DHA are shown in green for the closed conformation, and in orange and red, respectively, for the forred, respectively, for the open conformationtheof open conformation the MD replica 1. of the MD replica 1. open conformation of the MD replica 1. For the U-shaped DHA binding mode, the residues at the cavity bottom that mainly interact with the substrate are Phe352, Ala417, Cys559, and Gln590. For the more twisted DHA orientation in replica 1, the bottom cavity is formed by the sidechains of Ile357, His425, and Ile408. Those residues at the cavity bottom are shown in Figure 8. The changes Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 10 of 32 Int. J. Mol. Sci. 2023, 24, 6064 9 of 29 For the U-shaped DHA binding mode, the residues at the cavity bottom that mainly interact with the substrate are Phe352, Ala417, Cys559, and Gln590. For the more twisted in the DHA orientation in replica 1, the bottom cavity is formed by the sidechains of Ile357, interactions between the substrate’s tail and the protein residues at the cavity bottom His425, and Ile408. Those residues at the cavity bottom are shown in Figure 8. The changes couldinexplain the interactionshALOX12 why between themight onlytailpartially substrate’s follow and the protein the triad residues at thehypothesis. cavity bot- However, kinetictomexperiments could explain with A417I and why hALOX12 V418M might mutants only partially would follow behypothesis. the triad necessary to confirm this How- prediction that has ever, kinetic alreadywith experiments been shown A417I for AA. and V418M Finally, mutants it isbeworth would mentioning necessary to confirm that Leu407 is alsothis prediction a critical that has residue already been defining the shown for AA. U-shaped Finally, mode binding it is worth mentioning of DHA that in hALOX12, as has Leu407 is also a critical residue defining the U-shaped binding mode of DHA in hALOX12, been described for AA [22]. as has been described for AA [22]. 8. Main FigureFigure interactions 8. Main interactions between between thethe terminal terminal methyl methyl of DHAofandDHA andofresidues residues ofofthe bottom of the bottom hALOX12 0 s cavity hALOX12′s cavity forfor theclose the close (upper, (upper, blue) andand blue) openopen (lower, orange)orange) (lower, conformations. conformations. Here, we also analyzed the evolution of the hydrogen atoms H9 , H12 , and H15 , which are candidates to be abstracted in the first step of the catalytic mechanism. That is, we recorded the distances from the oxygen atom in the Fe(III)-OH− cofactor to the closest H9 , H12 , and H15 atoms (attached to C9 , C12 , and C15 , respectively) along the MD trajectory (see Figure 9). During the MD replica 1, hydrogen H12 remains at precatalytic distances (that is, smaller than 4 Å), although with some structural fluctuations. H9 also remains quite stable at precatalytic distances from the cofactor but moves farther away when the PLAT domain separates. The average H12 -OH− and H9 -OH− distances before the PLAT movement are 3.60 Å and 3.16 Å, respectively. The smaller average H9 -OH− distance compared to the H12 -OH− one does not agree with the regioselectivity of hALOX12, which gives 14-H(p)DHA as the main product. However, we have shown in previous mechanistic studies on the reactivity of other ALOXs that interatomic distances alone cannot explain H12, and H15 atoms (attached to C9, C12, and C15, respectively) along the MD trajectory (see Figure 9). During the MD replica 1, hydrogen H12 remains at precatalytic distances (that is, smaller than 4 Å), although with some structural fluctuations. H9 also remains quite stable at precatalytic distances from the cofactor but moves farther away when the PLAT domain separates. The average H12-OH− and H9-OH− distances before the PLAT movement Int. J. Mol. Sci. 2023, 24, 6064 10 of 29 are 3.60 Å and 3.16 Å, respectively. The smaller average H9-OH− distance compared to the H12-OH− one does not agree with the regioselectivity of hALOX12, which gives 14- H(p)DHA as the main product. However, we have shown in previous mechanistic studies on the of the molecular origin reactivity enzyme of regiospecificity. other ALOXs that interatomic As can be distances observedalone cannot9,explain in Figure H15 isthe mo- lecular origin of enzyme regiospecificity. As can be observed more distant from the cofactor than H12 and H9 when the two-domain protein is closed in Figure 9, H 15 is more dis- but tant from the cofactor than H 12 and H9 when the two-domain protein is closed but ap- approaches when the two domains open. The evolution of those three distances indicates proaches when the two domains open. The evolution of those three distances indicates that the opening of the two-domain protein is correlated with a movement of DHA to the that the opening of the two-domain protein is correlated with a movement of DHA to the entrance of the cavity, as mentioned above. Along the MD replica 2, the three distances entrance of the cavity, as mentioned above. Along the MD replica 2, the three distances maintain more stable maintain values morein accordance stable values inwith the stability accordance with the ofstability the DHA binding of the mode inmode in DHA binding this trajectory. In this case, theInaverage − this trajectory. this case,Hthe 12 -OH averagedistance H12-OHis the smallest − distance is the (2.99 Å) (2.99 smallest followed Å) followed by the average by H9the-OH − distance (3.44 Å), and − average H9-OH− distance (3.44the Å),H 15 -OH and the H15one -OH−(4.49 Å), in one (4.49 Å),agreement in agreement with with the regiospecificity observedobserved the regiospecificity for hALOX12. for hALOX12. Figure 9. Evolution with time of the Htime − , H -OH− , and H -OH− distances along the MD Figure 9. Evolution with 12 -OH of the H912-OH−, H9-OH−, 15 and H15-OH− distances along the MD replica replica 1 and replica 2 in 1 and hALOX12. replica 2 in hALOX12. As mentioned above, ALOXsabove, As mentioned are also enzymes ALOXs that are also are highly enzymes stereospecific. that are ALOX12 highly stereospecific. ALOX12 catalyzes the formation of 14S-hydroperoxy-4Z,7Z,10Z,12E,16Z,19Z-docosahexaenoic catalyzes the formation of 14S-hydroperoxy-4Z,7Z,10Z,12E,16Z,19Z-docosahexaenoic acid (14S-H(p)DHA),acid which (14S-H(p)DHA), which is the intermediate that is the13S,14S-epoxy-4Z,7Z,9E,11E,16Z,19Z- forms intermediate that forms 13S,14S-epoxy- docosahexaenoic acid (13S,14S-epoxy-DHA). This epoxide is then enzymatically hydrolyzed to form the final maresin product (MaR1, 7R,14S-dihydroxy-4Z,8E, 10E,12Z,16Z,19Z- docosahexaenoic acid). The stereochemistry of the bonds around C12 in 14S-H(p)DHA is critical for the formation of MaR1. We analyzed the MD trajectories of the DHA/hALOX12 complex and concluded that configurations of DHA with the angle between planes Π(C10C11C12) and Π(C12 C13 C14 ) smaller than 90◦ are prompted to give the E stereochemistry of the ∆12 bond in 14S-H(p)DHA. Most of the configurations of the MD replica 1 before the PLAT movement accomplish that structural condition, and there is also a fraction of those con- figurations when the two domains are separated (see Figure 10). In contrast, many DHA configurations present angles between the Π(C7 C8 C9 ) and Π(C9 C10 C11 ) planes with val- ues above 90◦ . These configurations would lead to a Z stereochemistry of the ∆9 bond in 11S-H(p)DHA. between planes Π(C10C11C12) and Π(C12C13C14) smaller than 90° are prompted to give the E stereochemistry of the Δ12 bond in 14S-H(p)DHA. Most of the configurations of the MD replica 1 before the PLAT movement accomplish that structural condition, and there is also a fraction of those configurations when the two domains are separated (see Figure 10). In contrast, many DHA configurations present angles between the Π(C7C8C9) and Int. J. Mol. Sci. 2023, 24, 6064 11 of 29 Π(C9C10C11) planes with values above 90°. These configurations would lead to a Z stereo- chemistry of the Δ9 bond in 11S-H(p)DHA. Figure 10. Evolution with10.time Figure of the with Evolution angletime between planesbetween of the angle formedplanes by C7formed C8 C9 and C79CC8C109 C by C 11 (in and C9Cred) 10C11 (in red) and by C C C and C C C (in blue) along the MD replica and by C10 C11 C12 and C12 C13 14 (in blue) along the MD replica 1 in hALOX12. 10 11 12 12 13 14 1 in hALOX12. 2.2. QM/MM Calculations 2.2. QM/MMof DHA Hydroperoxidation Calculations Catalyzed byCatalyzed of DHA Hydroperoxidation hALOX12 by hALOX12 2.2.1. Hydrogen2.2.1. Abstraction Reactions Hydrogen Catalyzed Abstraction Reactionsby hALOX12 Catalyzed by hALOX12 In Table 1 and Table S1, the In Tables QM/MM 1 and S1, the results QM/MM corresponding to the first results corresponding tostep of the the first stepcatalytic of the catalytic mechanism are givenmechanism are given (see Figure 11). (see We Figure 11). We have have calculated thecalculated potential the potential energy energy profiles profiles for for the H12, H9, and H15the H12, H9, and abstraction H15 abstraction reactions using asreactions using as initial initial structures structures several selected several snapshotsselected of snap- the MD trajectoryshots of the MD for replica trajectory 1. Those foror initial replica 1. Thosestructures precatalytic initial or precatalytic were filtered structures according were to filtered according to the following criteria. We searched precatalytic snapshots the following criteria. We searched precatalytic snapshots in which the initial distance between in which the initial distance between the three hydrogen atoms, H12, H9, and H15, and the oxygen atom of the the three hydrogen atoms, H12, H9, and H15, and the oxygen atom of− the hydrogen acceptor, hydrogen acceptor, that is, the oxygen of the Fe(III)-OH cofactor, was smaller than 4.0 Å. that is, the oxygen of the Fe(III)-OH− cofactor, was smaller than 4.0 Å. Table 1. Number of the MD frame used as initial coordinates for the QM/MM optimization in hALOX12, optimized H12proS -OH− and H9proR -OH− distances (in Å) at reactants, potential energy barriers, and reaction energies (in kcal/mol) for the H12proS and H9proR abstraction processes. The stereochemistry of the pentadienyl radicals centered at C12 and C9 is also given. Exponential average (in kcal/mol) energy barriers [46] for the two H-abstractions are included in the last row. H12proS H9proR a Pentadienyl b Pentadienyl Frame dreact H−OH− ∆E‡ ∆E dreact ∆E‡ ∆E Stereochemistry H−OH− Stereochemistry 322 3.5 15.9 −18.3 ZE 3.5 36.4 −13.9 ZZ 2253 5.4 20.5 −19.0 ZE 3.5 23.4 −19.1 ZZ 4986 5.4 19.3 −17.0 ZE 5.0 35.0 −14.1 ZZ 6993 3.9 18.0 −18.5 ZE 3.1 29.7 −14.5 ZZ 8721 5.0 21.9 −17.3 ZE 3.3 22.7 −13.6 ZZ 10,106 5.5 24.3 −16.8 ZE 3.4 25.6 −17.4 ZZ 12,860 5.2 22.9 −18.5 ZE 3.9 29.2 −15.3 ZZ 14,423 4.1 21.8 −19.1 ZE 3.4 27.5 −13.6 ZZ 18,168 3.7 17.1 −17.0 ZE 3.3 31.8 −13.8 ZZ 19,729 3.4 17.1 −17.0 ZE 3.9 43.7 −18.4 ZZ ‡ ∆EAV 17.1 23.9 a b H stands for H12proS ; H stands for H9proR . Int. J. Mol. Int. Sci. Sci. 24, 2023, J. Mol. 6064 2023, 24, x FOR PEER REVIEW 13 of 3212 of 29 Figure 11. Scheme of the DHA hydroperoxidation mechanism. PCET stands for proton coupled Figure 11. Scheme of the DHA hydroperoxidation mechanism. PCET stands for proton coupled electron transfer. electron transfer. Then, several of those precatalytic snapshots were optimized according to the Then, several QM/MM model of those precatalytic described snapshotssection, in the methodology were optimized according to minima and the corresponding the QM/MM model weredescribed located. Fromin the methodology those minima, thesection,potential and the corresponding energy minima were profile for the H-abstraction was located. From those minima, calculated using the the potentialbetween difference energythe profile CX-Hfor the H-abstraction X distances was calculated and the corresponding HX- using theOHdifference − ones as between the reaction thecoordinate. CX -HX distances In Table 1,andthe the corresponding results HX -OH−abstrac- for the most favorable ones as the reaction coordinate. tion reactions, the HIn12proS Table and 1, Hthe results 9proR for theare abstractions, most favorable collected. abstraction In Table S1, we havereactions, in- the cluded H12proS andtheH results 9proR for the rest abstractions, of the are H-abstraction collected. energy In Table profiles. S1, we The have values of included the H the proX - results for OH distances the− rest of the in H-abstraction Table 1 are givenenergy for the optimized profiles. The reactants valuesof the of H the12 and HproXH9 -OH − distances abstraction processes. in Table 1 areAsgiven can beforobserved, the QM/MM the optimized optimized reactants of thedistances H12 and are H longer than the initial 9 abstraction processes. As can be observed, the QM/MM optimized distances are longer than the there distances of the precatalytic structures. As for the potential energy barriers, initialisdistances no correlation between their values and those of the optimized H12proS-OH− or H9proR-OH− dis- of the precatalytic structures. As for the potential energy barriers, there is no correlation tances. Thus, the H9proR-OH− distances are shorter in almost all the optimized reactants, between their values and those of the optimized H12proS -OH− or H9proR -OH− distances. but the individual barriers− distancesare higher for the H9-abstraction than for the H12-abstraction. Thus, the H9proR -OH At the optimized are shorter in almost all the optimized reactants, but the radical products of both H-abstractions, the corresponding pentadienyl individual group is planar, and the electronic the barriers are higher for H9 -abstraction density is delocalizedthan over for fivethe carbonH12atoms. -abstraction. In DHA,At the optimized radical products of both H-abstractions, the corresponding the set of carbon atoms around C12 that will form the pentadienyl radical after H12-abstrac- pentadienyl group is planar, tion is closer to planarity than the set of carbon atoms around C9, which will form thethe set and the electronic density is delocalized over five carbon atoms. In DHA, of carbon pentadienylatoms around radical afterCH that will form 129-abstraction. Thethe pentadienyl larger radical the geometrical afterfrom change H12 -abstraction an ini- is closer to planarity tially nonplanar than the structure to aset of carbon planar atomsand, pentadienyl around C9 , which especially, will form the greater steric the hin-penta- drance dienyl to thatafter radical motion, the higher the contribution H9 -abstraction. The larger to thethegeometrical potential energy changebarrierfrom[24].an Bothinitially abstraction nonplanar processes structure toare exoergic, a planar as has already pentadienyl been and, obtained in especially, other the ALOXs greater withhindrance steric AA. All inmotion, to that all, the exponential the higher average energy barrier the contribution is 6.8 to the kcal/molenergy potential lower for H12-abstraction barrier [24]. Both ab- than for H 9-abstraction. This result agrees with the regioselectivity of hALOX12, which straction processes are exoergic, as has already been obtained in other ALOXs with AA. favors the abstraction of H12, leading to the formation of 14S-H(p)DHA after the addition All in all, the exponential average energy barrier is 6.8 kcal/mol lower for H12 -abstraction of O2 at C14 (see next section). In Table 1, the pentadienyl stereochemistry is given at both than for H9 -abstraction. This result agrees with the regioselectivity of hALOX12, which product radicals. It is remarkable that the stereochemistry of all the pentadienyl radicals favors fromtheC10abstraction to C14 (ZE) isof theHone 12 , leading needed for to the formation leading of 14S-H(p)DHA to the fully after the addition characterized 14S-hydroper- of O 2 at C14 (see next section). In Table 1, the oxo-4Z,7Z,10Z,12E,16Z,19Z-docosahexaenoic acid pentadienyl (see next section). stereochemistry (ZE) stands for is the given at both product radicals. stereochemistry in theItpentadienyl is remarkable thatthat radical the will stereochemistry become 10Z, of andall12Ethe inpentadienyl 14S- radicals H(p)DHA. fromHence, C10 tothe C14 (ZE) isofthe geometry theone QM/MMneeded for leading optimized products to the fully ratifies thecharacterized analysis 14S-hydroperoxo-4Z,7Z,10Z,12E,16Z,19Z-docosahexaenoic made of the most stable configurations of DHA along the MD trajectory acid (see(See nextFigure section). 10). (ZE) In contrast, stands for thethe stereochemistryinofthe stereochemistry thepentadienyl pentadienyl radicalradicalcentered that will at C 9 is (ZZ). become In this 10Z, and 12E case, the stereochemistry in 14S-H(p)DHA. Hence, of the11-H(p)DHA, geometry of being thethe minor product, QM/MM optimized has not been reported. products ratifies the In Figure analysis made 12, the of structures the most of the optimized stable reactant,oftransition configurations DHA along state structure, the MD and prod- (See trajectory uct of the H 12proS-abstraction initiated from snapshot 8721 are depicted. Figure 10). In contrast, the stereochemistry of the pentadienyl radical centered at C is (ZZ). 9 In this case, the stereochemistry of 11-H(p)DHA, being the minor product, has not been reported. In Figure 12, the structures of the optimized reactant, transition state structure, and product of the H12proS -abstraction initiated from snapshot 8721 are depicted. Int. J.Int. Mol. Sci. Sci. J. Mol. 2023, 24,24, 2023, 6064 x FOR PEER REVIEW 15 of 32 13 of 29 Figure Figure 12. Imagesofofthe 12. Images theoptimized optimized reactant reactant (a), transition (a), transition state structure state structure (b), and(b), and product product (c) of the (c) of the HH12proS 12proS abstraction abstraction initiated from initiated snapshot from 8721 in snapshot hALOX12. 8721 The distance in hALOX12. between Hbetween The distance 12proS and Hthe 12proS and oxygen atom in OH − is shown for the optimized reactant (a). At the transition state structure, the − the oxygen atom in OH is shown for the optimized reactant (a). At the transition state structure, distances between the shifting hydrogen with respect to C12 (donor atom) and the oxygen atom in the OHdistances − (acceptorbetween atom) arethe shifting depicted (b).hydrogen with The distance respect between thetoHCatom 12 (donor in the atom) nascentand waterthemole- oxygen atom incule OH − and (acceptor atom) C12 is plotted at theare depicted optimized (b). The product (c). distance between the H atom in the nascent water molecule and C12 is plotted at the optimized product (c). On the other hand, the H15 -abstraction (Table S1) has a low probability of occurring since the percentage of precatalytic structures (10%, see Figure 9) is quite lower than for the H12 (74%) and H9 (73%) abstractions. Int. J. Mol. Sci. 2023, 24, 6064 14 of 29 2.2.2. O2 Addition and Retro-Hydrogen Abstraction Reactions Catalyzed by hALOX12 The second step of the overall DHA hydroperoxidation process that leads to the final products consists of adding an oxygen molecule to the C10 –C14 or C7 –C11 DHA pentadienyl radicals, formed once H12 or H9 has been abstracted, respectively. Here, we only calculated the reaction pathway for O2 addition to the delocalized C10 –C14 radical that gives the main hydroperoxide product, 14S-H(p)DHA. We have explored different initial locations for the O2 molecule around C14 at the conformation of the DHA pentadienyl radical obtained from the different MD snapshots. The O2 molecules around C14 , taken as the origin of coordinates, have been initially placed at 3.0 Å along the x, y, and z Cartesian axes and along the bisector axes contained in the xy, xz, and yz planes. The structures selected were chosen by a visual analysis. This means that all the O2 molecules with close contacts or clashes with other residues were discarded. Then, QM/MM single-point energy calculations were carried out for the selected positions, and the higher energy structures were discarded. The most stable structures were optimized and then taken as starting points to build the reaction path for the oxygen addition to C14 . The addition pathway was calculated in the forward and backward directions using the distance from the attacking oxygen of the O2 molecule to C14 as the reaction coordinate. From the last structure of the backward path, the addition reactant structure was optimized. The O-C14 distance values for the converged QM/MM minimum energy geometries are presented in Table 2. All the oxygen addition reaction pathways present low potential energy barriers. The chirality of the oxygenated product and the geometry of the addition pathway are also included in Table 2. We have obtained three converged addition pathways corresponding to antarafacial additions leading to peroxyl radicals with S stereochemistry at C14 , in agreement with the experimental stereochemistry assigned to the final hydroperoxide product, 14S-H(p)DHA. In Figure 13, the reactant, transition state structure, and the product of the addition process to C14 , initiated from frame 8721, are depicted. The final step of the overall hydroperoxidation mechanism consists of a retro-hydrogen abstraction from the Fe(II)-H2 O cofactor to the peroxyl DHA radical. This reaction leads to the final 14S-H(p)DHA product, and hALOX12 recovers the initial Fe(III)-OH− cofactor’s state. However, as the oxygen molecule has been added following an antarafacial approach with respect to the cofactor, the carbon chain of the peroxyl radical needs to rotate before the proper retro-hydrogen transfer can take place. We calculated the potential energy profiles for this carbon chain rotation using the C13 C14 C15 C16 dihedral as the reaction coordinate for the three pathways initiated at snapshots 8721, 10,106, and 18,168. The potential energy barriers measured from the corresponding previous minima (that is, the addition products) on the QM/MM potential energy surface for the three rotations range from 5.8 to 6.6 kcal/mol. In Figure 14, the reactant, transition state structure, and the product of the rotation process, initiated from frame 8721, are presented. Note that the reactant minimum in Figure 14a has the same structure as the addition product minimum in Figure 13c but is plotted from a different perspective to better follow the carbon chain rotation. Table 2. Number of the MD frame used as initial coordinates for the QM/MM optimization in hALOX12; optimized O-C14 distances (in Å) at the reactants and potential energy barriers for the addition reaction as well as the chirality of the peroxyl product and the geometry of the addition approach; potential energy barriers for the rotation of the peroxyl radical’s carbon chain and for the retro-hydrogen abstraction (reorganization of the peroxyl radical and retro-hydrogen abstraction itself). All energies are in kcal/mol. ‡ Chirality of Geometry of ‡ ‡ ‡ Frame dreact O−C14 ∆EO ∆ECrot ∆Ereorg ∆Eretro−Habs 2 add Product Addition 8721 3.2 3.0 S Antarafacial 6.6 10.2 33.3 10,106 3.2 4.3 S Antarafacial 5.8 - - 18,168 3.1 3.1 S Antarafacial 5.9 13.1 20.6 Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 17 of 3 Int. J. Mol. Sci. 2023, 24, 6064 15 of 29 Figure13. Figure 13.Images Images of of thethe optimized optimized reactant reactant (a), transition (a), transition state structure state structure (b), and(b), and (c) product product of the (c) of th OO22addition additionatat C14 C14 (frame (frame 8721) 8721) in hALOX12. in hALOX12. The distances The distances between between C14 andC 14 and the the attacking attacking oxygen oxygen atomare atom aregiven given forfor thethe reactant reactant andand transition transition state structures. state structures. The final step of the overall hydroperoxidation mechanism consists of a retro-hydro gen abstraction from the Fe(II)-H2O cofactor to the peroxyl DHA radical. This reaction leads to the final 14S-H(p)DHA product, and hALOX12 recovers the initial Fe(III)-OH cofactor’s state. However, as the oxygen molecule has been added following an anta rafacial approach with respect to the cofactor, the carbon chain of the peroxyl radica needs to rotate before the proper retro-hydrogen transfer can take place. We calculated the potential energy profiles for this carbon chain rotation using the C13C14C15C16 dihedra (that is, the addition products) on the QM/MM potential energy surface for the three ro tations range from 5.8 to 6.6 kcal/mol. In Figure 14, the reactant, transition state structure and the product of the rotation process, initiated from frame 8721, are presented. Note that the reactant minimum in Figure 14a has the same structure as the addition produc Int. J. Mol. Sci. 2023, 24, 6064 minimum in Figure 13c but is plotted from a different perspective to better 16 of follow 29 the carbon chain rotation. Figure14. Figure 14.Images Images of of thethe optimized optimized reactant reactant (a), transition (a), transition state structure state structure (b), and(b), and (c) product product of the (c) of th carbonchain carbon chainrotation rotation of the of the peroxyl peroxyl radical radical (frame (frame 8721) 8721) in hALOX12. in hALOX12. From the rotated peroxyl radical the retro-hydrogen abstraction from the Fe(II)-H2 O cofactor to the peroxyl DHA radical can be initiated. We calculated the corresponding QM/MM potential energy profile versus a reaction coordinate defined as the difference between the breaking H-OH (water ligand) bond distance and the forming H-OOC14 one. This abstraction process involves two different structural steps. First, there is a reorgani- From the rotated peroxyl radical the retro-hydrogen abstraction from the Fe(II)-H2O cofactor to the peroxyl DHA radical can be initiated. We calculated the correspondin QM/MM potential energy profile versus a reaction coordinate defined as the differenc Int. J. Mol. Sci. 2023, 24, 6064 between the breaking H-OH (water ligand) bond distance and the forming17H-OOC of 29 14 one This abstraction process involves two different structural steps. First, there is a reorgan zation of the peroxyl radical that correlates with the rotation of the peroxo group at C from of zation thethe antarafacial to the peroxyl radical suprafacial that side.the correlates with Therotation potential energy of the peroxo barriers group at given C14 in Tabl 2 for this structural reorganization (10.2 and 13.1 kcal/mol for the pathways from the antarafacial to the suprafacial side. The potential energy barriers given in Table 2 initiated from snapshots for 8721 reorganization this structural and 18,168, respectively) (10.2 and 13.1are calculated kcal/mol for thefrom the corresponding pathways initiated from reactan snapshots minima, which are the products obtained from the rotation pathway. Thereactant 8721 and 18,168, respectively) are calculated from the corresponding reorganizatio minima, which aretothe corresponding products10,106 snapshot obtaineddidfrom notthe rotation pathway. converge. The transitionThe reorganization state structure (TS1 corresponding and the product to snapshot geometry10,106 ordid not converge. intermediate The transition (INT1) state structure (TS1) of this reorganization and are de process the product geometry or intermediate (INT1) of this reorganization process are depicted in picted in Figure 15 for the reaction initiated from snapshot 8721. The final step of th Figure 15 for the reaction initiated from snapshot 8721. The final step of the mechanism mechanismtocorresponds corresponds the movement toofthe themovement of the shifting shifting hydrogen. hydrogen. The abstraction The energy potential abstraction po tential energy barriers calculated from the optimized intermediate barriers calculated from the optimized intermediate are 33.3 and 20.6 kcal/mol for are 33.3 the and 20. kcal/molinitiated pathways for the from pathways framesinitiated 8721 and from 18,168,frames 8721 and respectively. 18,168, In Figure 15, respectively. the plots of theIn Figur 15, the plots ofstructures stationary-point the stationary-point structures (transition state (TS2) and (transition product) of state this(TS2) and are final step product) also of thi depicted for the pathway of snapshot 8721. final step are also depicted for the pathway of snapshot 8721. Figure 15. Cont. Int.Int. J. Mol. Sci. J. Mol. Sci.2023, 2023,24, 24,x6064 FOR PEER REVIEW 18 of 29 20 of 3 Figure15. Figure 15.Images Imagesof of thethe transition transition statestate structure structure (TS1) (TS1) (a), intermediate (a), intermediate (INT1) (INT1) (b), (b), transition transition state stat structure (TS2) (c), and the final 14S-H(p)DHA product (d) (frame 8721) in hALOX12. structure (TS2) (c), and the final 14S-H(p)DHA product (d) (frame 8721) in hALOX12. The distances The distance betweenthe between the outer outer oxygen oxygen of the of the peroxo peroxo groupgroup and and the the water water hydrogen hydrogen atom areatom givenare forgiven for TS1 an TS1 and INT1. The distances corresponding to the shifting hydrogen atom are indicated for TS2 and the final the fina INT1. The distances corresponding to the shifting hydrogen atom are indicated for TS2 and 14S-H(p)DHA 14S-H(p)DHA product. product. InInFigure Figure 16,16, wewe have depicted have an overall depicted energy an overall scheme energy of the mechanistic scheme steps of the mechanistic step initiated from the reactant of the H -abstraction reaction (frame 8721) initiated from the reactant of the H12proS-abstraction reaction (frame 8721) bound to th 12proS bound to the solvated solvatedenzyme enzyme plus an an plus oxygen oxygenmolecule within molecule the water within box at box the water the entrance of the of th at the entrance enzyme cavity. The position of the oxygen molecule has been optimized at 12.1 Å from the enzyme cavity. The position of the oxygen molecule has been optimized at 12.1 Å from C14 atom of the substrate. This structure has been taken as the zero of energy of the overall the C14scheme. energy atom ofThe theapproach substrate. Thisoxygen of the structure has been molecule to 3.1taken Å from asCtheatzero of energy of th the addition 14 overall structure reactant energy scheme. representsThe approach a falloff of the in energy oxygen of 53.3 molecule kcal/mol to 3.1 to with respect Å the from C14 at H12proS - the ad abstraction product. Consequently, the barriers for the O2 addition, carbon chain rotation, to th dition reactant structure represents a falloff in energy of 53.3 kcal/mol with respect H12proS and -abstractionabstraction retro-hydrogen product. Consequently, are very far below thethebarriers referencefor the O2 addition, structure. carbon chain Therefore, none of those processes rotation, can be the rate-determining and retro-hydrogen abstraction are step veryoffarthe hydroperoxidation below reaction. There the reference structure. In contrast, fore, none the H12proSprocesses of those -abstraction cantransition state structure is 17.3 be the rate-determining stepkcal/mol above the of the hydroperoxidation reference structure (see Table 1 and Figure 16), so it is the rate-determining reaction. In contrast, the H12proS-abstraction transition state structure is 17.3 step. Here, we abov kcal/mol are assuming that the oxygen molecule approaches the enzyme once the H-abstraction the reference structure (see Table 1 and Figure 16), so it is the rate-determining step. Here process has taken place. we are assuming that the oxygen molecule approaches the enzyme once the H-abstraction process has taken place. Int. J. Mol. Int. Sci. Sci. 24, 2023, J. Mol. 6064 2023, 24, x FOR PEER REVIEW 21 of 3219 of 29 Figure Overall 16.16. Figure Overallenergy energyscheme scheme of the the HH12proS -abstraction, -abstraction, 12proS oxygen oxygen addition, addition, carbon carbon chain rotation, chain rotation, andand retro-hydrogenabstraction retro-hydrogen abstraction steps steps of of the theDHA DHAhydroperoxidation hydroperoxidationmechanism mechanismin hALOX12. All in hALOX12. All energies are in kcal/mol. The zero of energies corresponds to the reactant of the H12proS abstraction energies are in kcal/mol. The zero of energies corresponds to the reactant of the H12proS abstraction step bound to the solvated enzyme plus an oxygen molecule within the water box whose position stephas bound to the solvated been optimized at 12.1 enzyme plus Å from the an oxygen substrate’s C14. molecule within the water box whose position has been optimized at 12.1 Å from the substrate’s C14 . 2.3. MD Simulations of the DHA/pigALOX15-mini-LOX Complex 2.3. MD Simulations of the DHA/pigALOX15-mini-LOX Complex In this section, the results corresponding to the two MD trajectories of the solvated In this section, the resultscomplex DHA/pigALOX15-mini-LOX corresponding to thefor are presented twotheMD trajectories sake of thewith of comparison solvated DHA/pigALOX15-mini-LOX DHA/hALOX12 behavior. In Figure complex S4, weare presented have plotted thefor the sake protein of comparison backbone RMSD ver- with DHA/hALOX12 behavior. sus time for the two In Figure MD replicas. S4, we The results have show plotted a stable the protein conformation backbone of the C-terminalRMSD versus domaintime for the two MD replicas. of pigALOX15-mini-LOX, The results even though show adomain the N-terminal stableisconformation not present. Theof the carbon-chain C-terminal domainRMSDof of DHA also indicates a rather pigALOX15-mini-LOX, evensteady though binding mode alongdomain the N-terminal the two is not trajectories. present. The carbon-chain RMSD of DHA also indicates a rather steady binding mode along theIn the twopigALOX15-mini-LOX trajectories. active site, DHA also adopts a U-shaped binding mode like that described in hALOX12. However, the carboxylate head of DHA does not interact In the pigALOX15-mini-LOX active site, DHA also adopts a U-shaped binding mode like with Arg403 (Arg402 in hALOX12), as can be observed in Figure 17 along the two MD that described in hALOX12. However, the carboxylate head of DHA does not interact with trajectories. Arg403 (Arg402 in hALOX12), as can be observed in Figure 17 along the two MD trajectories. This interaction is lost because a displacement of the α2-helix in pigALOX15-mini-LOX moves away from the Glu176 . . . Arg403 pair (Glu175 . . . Arg402 in hALOX12) from the cavity entrance. Moreover, pigALOX15-mini-LOX has a Gly in position 407 instead of the Gln406 of hALOX12, which does not interact with Arg403 and cannot maintain this arginine close to the DHA carboxylate. Therefore, in its location in pigALOX15-mini-LOX, the carboxylate head of DHA only interacts with some water molecules. In addition, DHA also establishes π–π interactions between its double bonds and aromatic residues of the active site that stabilize the substrate binding mode. The stacking interactions between Phe175 and ∆4 and ∆7 are plotted in Figure 18. Along the MD replicas, Phe175 exchanges its interaction between those two double bonds, except during the first 30 ns of replica 2 Int. J. Mol. Sci. 2023, 24, 6064 20 of 29 Int. J. Mol. Sci. 2023, when 24, xthe FORstacking PEER REVIEWwith both double bonds at the same time. The rest of the22π–π is present of 32 interactions presented in Figure S5 are weaker in pigALOX15-mini-LOX than in hALOX12. Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 23 of 32 Figure 17. Distances of17. Figure a carboxylate’s oxygen atom Distances of a carboxylate’s of DHA oxygen toDHA atom of the closest hydrogen to the closest atom hydrogen in Arg403 atom in Arg403 versus versus time for the MDtime for the replica 1 MD andreplica 1 and replica 2 inreplica 2 in pigALOX15-mini-LOX. pigALOX15-mini-LOX. This interaction is lost because a displacement of the α2-helix in pigALOX15-mini- LOX moves away from the Glu176...Arg403 pair (Glu175…Arg402 in hALOX12) from the cavity entrance. Moreover, pigALOX15-mini-LOX has a Gly in position 407 instead of the Gln406 of hALOX12, which does not interact with Arg403 and cannot maintain this argi- nine close to the DHA carboxylate. Therefore, in its location in pigALOX15-mini-LOX, the carboxylate head of DHA only interacts with some water molecules. In addition, DHA also establishes π–π interactions between its double bonds and aromatic residues of the active site that stabilize the substrate binding mode. The stacking interactions between Phe175 and Δ4 and Δ7 are plotted in Figure 18. Along the MD replicas, Phe175 exchanges its interaction between those two double bonds, except during the first 30 ns of replica 2 when the stacking is present with both double bonds at the same time. The rest of the π– π interactions presented in Figure S5 are weaker in pigALOX15-mini-LOX than in hA- LOX12. Figure 18. Distances between Figure Phe175 18. Distances sidechain between Phe175and the closest sidechain and theatom the ∆ ofatom closest 4 and4∆7 double of the bonds Δ and Δ7 double bonds versus time for the MD replica 1 and for replica 2 in pigALOX15-mini-LOX. versus time for the MD replica 1 and for replica 2 in pigALOX15-mini-LOX. As for H12-OH−, H9-OH−, and H15-OH− distances along the MD replica 1 and replica 2, we have plotted them in Figure 19. H12 remains the closest to the OH- group of the cofactor the main part of the trajectory. H9 and H15 present a similar number of precatalytic structures when considering both MD replicas. However, when H9 is closer, H15 is farther Int. J. Mol. Sci. 2023, 24, 6064 21 of 29 As for H12 -OH− , H9 -OH− , and H15 -OH− distances along the MD replica 1 and replica 2, we have plotted them in Figure 19. H12 remains the closest to the OH- group of the cofactor the main part of the trajectory. H9 and H15 present a similar number of precatalytic structures when considering both MD replicas. However, when H9 is closer, H15 is farther away, and vice versa. Our MD simulation indicates that the DHA binding mode in pigALOX15-mini-LOX could lead to 14-H(p)DHA as a major product but also, with less probability, to 11-H(p)DHA and 17-H(p)DHA as minor products. Experimentally, Kühn Int. J. Mol. Sci. 2023, andPEER 24, x FOR coworkers REVIEW obtained a relative share of 76.4% for 14-H(p)DHA and 24 23.6% of 32 for 17-H(p)DHA for the Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW specificity of pigALOX15 [21]. 24 of 32 Figure 19. Evolution with time with of the Hof12the -OH − , H− -OH−− , and H -OH − distances along the MD Figure 19. Evolution time H12-OH ,9H9-OH , and H15-OH 15 − distances along the MD replica 1 and replica 2 in pigALOX15-mini-LOX. replica 1 and replica 2 in pigALOX15-mini-LOX. − Figure 19. Evolution with time of the H12-OH , H9-OH , and H15-OH distances along the MD replica − − 1 andWhen replicawe 2 in pigALOX15-mini-LOX. calculated the angle between planes Π(C10C11C12) and Π(C12C13C14) of DHA, When wethecalculated the angle between planes Π(C10 Figure values obtained are above 90° most of the time (see C11 C1220). ) and Π(C12 Cbefore, As indicated 13 C14 ) of DHA, the values obtained When we are above calculated the 90 ◦ most angle betweenof the planestime Π(C (see 10C11C12)Figure and Π(C 20). 12C13CAs indicated 14) of DHA, this result means that the configurations of DHA would give a Z stereochemistry of the the before, this result Δ bond 12 values means obtained are above 90° that the configurations in 14S-H(p)DHA. most of The angle between the of DHAtime (see planeswouldFigure Π(C7C8C give 20). As a Π(C 9) and indicated before, Z stereochemistry 9C10C11) of DHA of this result means all that the configurations of DHA would give a Z stereochemistry of the the ∆12 bond inbehaves similarly 14S-H(p)DHA. along The the trajectory. angle between planes Π(C 7 C 8 C 9 ) and Δ bond in 14S-H(p)DHA. The angle between planes Π(C7C8C9) and Π(C9C10C11) of DHA 12 Π(C 9 C 10 C 11 ) of DHA behaves behaves similarly all along similarly thethe all along trajectory. trajectory. Figure 20. Evolution with time of the angle between planes formed by C7C8C9 and C9C10C11 (in red) and by C10C11C12 and C12C13C14 (in blue) along the MD replica 1 in pigALOX15-mini-LOX. Figure 20. Evolution Figurewith time ofwith 20. Evolution the time angle between of the planesplanes angle between formed by Cby7 C formed C8C C87C 9 9and C99CC1010 and C C11C(in (in red) 11 red) and by C10 C11 Cand 12 by C10C11C12 and C12C13C14 (in blue) along the MD replica 1 in pigALOX15-mini-LOX. and C C C 12 13 14 (in blue) along the MD replica 1 in pigALOX15-mini-LOX. Int. J. Mol. Sci. 2023, 24, 6064 22 of 29 2.4. QM/MM Calculations of H12 -Abstraction from DHA Catalyzed by pigALOX15-mini-LOX In this section, we present the results of the QM/MM calculations corresponding to the H12 -abstraction. We have calculated the potential energy profiles for the H12 transfer using as initial structures several selected snapshots of the MD trajectory for replica 1. Following the methodology explained in the corresponding section, we calculated the potential energy profile of this catalytic step using, as a reaction coordinate, the difference between the C12 -H12 and the H12 -OH− distances. In Table 3, the potential energy barriers are given along with their exponential average. The obtained value of 17.6 kcal/mol is nearly the same as for the H12proS -abstraction in hALOX12. This theoretical result verifies, in agreement with experiments, that pigALOX15-mini-LOX is active and can catalyze the conversion of DHA to 14-H(p)DHA beginning with the abstraction of H12proS [21]. In Table 3 the stereochemistry of the pentadienyl radical centered at C12 is given. All the radical products present a (ZZ) stereochemistry, as the configurations of the MD trajectory had predicted. This is not such a stable stereochemistry for those pentadienyl radicals in the enzyme’s binding pocket as the (ZE) configuration obtained in hALOX12. Hence, the corresponding reaction energies are not as negative as in hALOX12 (Table 1). Moreover, as commented above, the (ZZ) stereochemistry is not the adequate configuration for the formation of 13S,14S-epoxy-4Z,7Z,9E,11E,16Z,19Z-docosahexaenoic acid and for finally leading to 7R,14S-dihydroxy-4Z,8E, 10E,12Z,16Z,19Z-docosahexaenoic acid (MaR1). For this reason, our hypothesis here is that pigALOX15-mini-LOX forms 14S-hydroperoxo- 4Z,7Z,10Z,12Z,16Z,19Z-docosahexaenoic acid, but this hydroperoxide would not give the bioactive SPM MaR1 that it is synthesized by hALOX12. Table 3. Optimized H12proS -OH− distances (in Å) at reactants, potential energy barriers, and reaction energies (in kcal/mol) for the H12proS abstraction process. The stereochemistry of the pentadienyl group around C12 at the product radical is also given. The exponential average potential energy barrier [46] (in kcal/mol) for the H12proS -abstraction is included in the last row. H12proS a Pentadienyl Frame dreact H−OH− ∆E‡ ∆E Stereochemistry 940 3.2 20.9 −12.5 ZZ 1548 3.0 19.8 −11.8 ZZ 1709 2.9 24.5 −12.4 ZZ 2132 2.8 22.1 −9.7 ZZ 2515 3.2 27.2 −8.9 ZZ 3020 2.8 19.8 −10.8 ZZ 3556 3.2 24.5 −9.9 ZZ 4502 2.7 15.9 −16.1 ZZ 5758 2.6 20.4 −12.1 ZZ 9207 2.6 22.6 −13.1 ZZ ‡ ∆EAV 17.6 a H stands for H12proS . 3. Materials and Methods 3.1. Protein Setup The model structure of human ALOX12 was taken from the AlphaFold Data Base (Unipro- tID P18054) [42]. For pigALOX15-mini-LOX, we used the x-ray structure launched at the Protein Data Bank with PDB code 3RDE [36], removing the ligand bound to the protein. The two structures were protonated through a web interface (www.playmolecule.org) (accessed Int. J. Mol. Sci. 2023, 24, 6064 23 of 29 on 19 January 2021) [47]. A pH = 7.0 was employed for the titratable residues. The protonation state for the iron coordination sphere was corrected by hand to ensure a correct description. 3.2. Molecular Docking Simulations Docking calculations were performed with GOLD5.8 [48] to obtain stable binding modes of DHA in the active site of hALOX12 and pigALOX15-mini-LOX. For the conforma- tional search of the ligand, the torsion angle distribution of the Cambridge Database was used. In contrast, some protein side chains were treated as flexible. The GOLD’s option to consider the interactions of organic ligands with metal ions in metalloenzymes was activated, but it restricted the docking exploration to hexacoordinated geometries of iron. The binding side was defined as a sphere of 20 Å around the iron atom. This conformational space was explored using the Lamarckian algorithm included in GOLD. One hundred DHA-binding poses were generated and grouped into clusters. We used the ChemScore fitness function to rank the docking poses and estimate their binding free energies. 3.3. Molecular Dynamics Simulations The best-ranked docking poses of DHA bound to hALOX12 and pigALOX15-mini- LOX were selected to initiate the MD simulations. For the solvated DHA/hALOX12 complex, we used the ff19SB force field [49] to calculate the potential energy of the protein and the OPC force field [50] to describe the water molecules, as recommended [49]. For the solvated DHA/pigALOX15-mini-LOX complex, we used the ff14SB force field [51] for the protein and the TIP3P force field to describe the water molecules as recommended [51]. We developed specific parameters for DHA using the AMBER standard protocol with Antechamber and Parmchk2 modules. The GAFF2 library [52] was used as the source for these parameters. The substrate structure was optimized employing the B3LYP/6-31G(d) level of theory. The atomic charges were set to fit the electrostatic potential generated at the B3LYP/6-31G(d) level of theory by the restrained electrostatic potential (RESP) model. The atomic charges were calculated according to the Merz–Kollman scheme [53] using Gaussian16 [54]. An unprotonated state was established for the DHA substrate in the two complexes, considering the physiological conditions. We also developed specific MM parameters for the iron atom and its first coordination sphere using the MCPB.py procedure [55] within the bonded model and the Seminario method for the force constant calculations [56]. For hALOX12, the iron ligands are His360, His365, His540, Asn544, Ile663, and the OH− group. For pigALOX15-mini-LOX, the iron ligands are His250, His255, His430, His434, Ile552, and the OH− group. The protocol recommended by the AMBER package using the tLeap module was used to assemble the DHA/hALOX12 and DHA/pigALOX15-mini-LOX systems, solvate those complexes with an orthorhombic box of pre-equilibrated water molecules with a buffer of 10 Å and neutralize the total charge by adding sodium cations. The resulting systems contain around 80,000 atoms, of which about 10,600 belong to the protein, for the DHA/hALOX12 system, and 64,000 atoms, of which about 8800 belong to the protein, for the DHA/pigALOX15- mini-LOX. The rest of the atoms correspond to water molecules and salt ions. The molecular dynamics (MD) simulations were run using either the AMBER20 (hALOX12) [57] or AMBER16 (pigALOX15-mini-LOX) [58] GPU (CUDA) version of the PMEMD package [59,60]. Initially, the systems were submitted to 110,000 energy mini- mization steps combining the steepest decent and conjugate gradient methods to remove close contacts. In the first 5000 steps, harmonic restraints were applied to the protein and substrate atoms with a force constant of 5.0 kcal mol−1 Å−2 so that only the solvent and ions were relaxed. In the following 5000 steps, harmonic restraints were applied to the protein backbone with the same force constant as before. Finally, the whole system was kept free during the last 100,000 steps. Then, MD simulations using periodic boundary conditions and the particle-mesh Ewald approach to introduce long-range electrostatic effects were performed. The DHA/hALOX12 system was gently heated using six 20 ps steps, incrementing the temperature by 50 K each step (0–300 K) under constant volume. Int. J. Mol. Sci. 2023, 24, 6064 24 of 29 For the DHA/pigALOX15-mini-LOX system, the temperature was increased by 10 steps of 30 K (0–300 K) during 20 ps each step. After heating, we calculated an MD trajectory of 1 ns within the NPT ensemble (300 K, 1 atm) to adjust the volume of the orthorhombic box and relax the density to a value of around 1 g cm−3 . During the heating and the compressing, harmonic restraints were applied to the protein backbone with a force con- stant of 5.0 kcal mol−1 Å−2 , whereas the rest of the system was kept free. The Langevin equilibration scheme [61] was used to control and equalize the temperature, while the pressure was adjusted by the Berendsen barostat [62]. Next, an equilibration stage of 10 ns at constant volume and temperature (300 K) was performed. Finally, we ran production MD trajectories of 250 ns and 100 ns for hALOX12 and pigALOX15-mini-LOX, respectively, within the same NVT ensemble. For each system, two MD replicas were calculated. A time step of 1 fs was used along the whole MD trajectories. Bonds involving hydrogen were constrained with the SHAKE algorithm [63]. The non-bonding interactions have been calculated with a cutoff of 9 Å. The analysis of the MD trajectories was carried out with a Python package and RCBS.py [64]. Matplotlib was used for plotting [65]. 3.4. QM/MM Calculations The modular program package ChemShell 3.7 [66,67] was employed to carry out the QM/MM computations. TURBOMOLE 7.0 [68] was used for the DFT calculations and the DL_POLY 5.0 [69] module in ChemShell for the MM part. The QM region was described by the B3LYP hybrid function [70]. The 6–31G(d) Pople basis set [71] was employed for the C, H, O, and N atoms, while the LANL2DZ basis set [72] was used for the Fe atom. The QM/MM partition is depicted in Figure 21. As for the hydrogen abstractions, the QM region was defined by all the atoms of the DHA substrate, which are found between C6 and C17 ; 11 atoms for each His residue of the iron coordination sphere (His360, His365, and His540 in hALOX12 and His250, His255 His430, and His434 in pigALOX15-mini-LOX); 8 atoms of the iron ligand Asn554 in hALOX12; and 3 atoms of the iron ligand Ile residue (Ile 663 in hALOX12 and Ile552 in pigALOX15-mini-LOX) and the FeIII -OH− cofactor. For oxygenations, the DHA carbon chain rotation, and the hydrogen retrodonation, this region was enlarged by an oxygen molecule. Seven link atoms were included to define the QM/MM boundary: five between the bonds Cα-QM atoms of the five residues in the iron coordination sphere and two bonded to the aliphatic carbon atoms of the lipid substrate (placed between C5 -C6 and C17 -C18 ). An electronic embedding scheme was employed to treat the interaction between the QM and MM subsystems. We also used the charge shift algorithm to minimize overpolarization effects. Cutoffs were not introduced to treat the nonbonding MM and QM/MM interactions [73]. The QM/MM optimizations have been carried out by employing the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm [74] for energy minimizations and scans of reaction pathways. For these optimizations to minima, a microiterative scheme [75] was considered using hybrid delocalized coordinates (HDLC) [76]. As for the transition- state searches, the dimer method was used [77]. These algorithms are implemented in the DL_FIND geometry optimization library [78] of ChemShell. For these QM/MM models, all water molecules outside a 17 Å radius volume centered on the ligand molecule were removed. The active region was defined by all residues and water molecules inside a 15 Å radius sphere centered on C12 of the ligand molecule. This region was allowed to move freely (≈2200 atoms), while the atoms left out were kept frozen during the optimization. Roughly 12,000 atoms were considered in the QM/MM calculations. The images of structures were plotted with UCSF CHIMERA [79]. Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 28 of 32 Int. J. Mol. Sci. 2023, 24, 6064 25 of 29 Figure 21. QM/MM partition in the DHA/hALOX12 (a) and DHA/pigALOX15-mini-LOX systems (b). Figure 21. QM/MM partition in the DHA/hALOX12 (a) and DHA/pigALOX15-mini-LOX systems (b). 4. Conclusions Human ALOX12 is a lipoxygenase that catalyzes the oxidation of the ω-3 fatty acid 4. Conclusions 4,7,10,13,16,19-docosahexaenoic acid (DHA), forming the hydroperoxides, 14S-H(p)DHA Human ALOX12asismajor and 11S-H(p)DHA, a lipoxygenase and minor that catalyzes products, the oxidation respectively. of the ω-3 The biological fatty acid relevance of 4,7,10,13,16,19-docosahexaenoic acid (DHA), forming the hydroperoxides, the corresponding reduced oxylipins, 14S-hydroxydocosahexaenoic acid (14S-HDHA) and 14S-H(p)DHA and 11S-H(p)DHA, as major and 11S-hydroxydocosahexaenoic acidminor products,has (11S-HDHA), respectively. The biological been associated relevance with the regulation of the corresponding reduced oxylipins, 14S-hydroxydocosahexaenoic of thrombus formation in vivo. PigALOX15 also catalyzes the conversion of DHA to acid (14S-HDHA) and 11S-hydroxydocosahexaenoic 14S-H(p)DHA, but the minor product aciddetected (11S-HDHA), has been associated with the regula- is 17S-H(p)DHA. tion of Inthrombus this paper,formation in vivo.out we have carried PigALOX15 docking and alsomolecular catalyzesdynamics the conversion of DHA simulations, to plus 14S-H(p)DHA, but the minor product detected is 17S-H(p)DHA. QM/MM calculations, to analyze the molecular details of the complete DHA hydroperoxi- dationIn mechanism this paper, we have carried in hALOX12. outthe For docking andthe first time, molecular protein dynamics simulations, setup of hALOX12 has plus been QM/MM calculations, to analyze the molecular details of the complete modeled using the structure provided by the AlphaFold server. For pigALOX15, the same DHA hydroperox- idation mechanism methodology in hALOX12. has been For the first used to calculate time, the the energy protein barrier of setup of hALOX12 the first step of thehas been catalytic modeled mechanism using thethe using structure providedof X-ray structure byan theN-terminal AlphaFoldtruncation server. Forvariant pigALOX15, the same of pigALOX15. methodology The DHAhas been used to peroxidation calculate reaction the energy is initiated barrier by the of the first of H-abstraction step H12,of H the catalytic 9 , or H15 of mechanism DHA by theusing the X-ray Fe (III)-OH − cofactor. structureNext, of anan N-terminal truncation oxygen molecule variant inserts on of pigALOX15. C14 , C11 , or C17 The DHA following peroxidation an antarafacial reactionFinally, approach. is initiated thereby is athe H-abstractiondonation retro-hydrogen of H12, Hfrom 9, or H 15 of the Fe DHA (II)-H2by the Fe (III)-OH O cofactor − to the peroxylcofactor. Next, radical, an oxygen leading molecule inserts to the corresponding on C14, C11, or hydroperoxides withC17S following an antarafacial stereochemistry. approach. The calculated QM/MM Finally, there energy is a retro-hydrogen barriers donation for each mechanistic stepfrom the confirm that Fe the reaction (II)-H process 2O cofactor to theis peroxyl viable inradical, the hALOX12leadingactive to thesite. Interestingly,hydroperoxides corresponding the stereochem- Int. J. Mol. Sci. 2023, 24, 6064 26 of 29 istry (ZE) of the pentadienyl radicals from C10 to C14 , formed after H12pros abstraction, is the one needed for leading to the fully characterized 14S-hydroperoxy-4Z,7Z,10Z,12E,16Z,19Z- docosahexaenoic acid because (ZE) in the pentadienyl radical becomes 10Z and 12E in 14S-H(p)DHA. Only with this stereochemistry can 14S-H(p)DHA be converted by ALOX12 to 7R,14S-dihydroxy-4Z,8E,10E,12Z,16Z,19Z-docosahexaenoic acid (MaR1). This fact is rele- vant because MaR1 is a potent SPM that actively participates in the inflammation-resolution phase of many diseases and has also been described as a novel antiplatelet agent. In pigALOX15-mini-LOX, our calculations support that 14S-H(p)DHA can be formed, but the stereochemistry (ZZ) of the pentadienyl radicals from C10 to C14 , obtained after H12pros abstraction, is not adequate for MaR1 formation. This is an excellent example of how an enzyme governs the stereochemistry of the catalytic mechanism that leads to a given bioactive product (MaR1 in this case) and how this stereochemistry control might be species specific. Supplementary Materials: The supporting information can be downloaded at: https://www.mdpi. com/article/10.3390/ijms24076064/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. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. 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. Steinhilber, D. Lipoxygenases: An Introduction. In Lipoxygenases in Inflammation; Steinhilber, D., Ed.; Springer: Cham, Switzerland, 2016; pp. 1–6; ISBN 978-3-319-27764-6. 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] 3. 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] 4. Hamberg, M.; Samuelsson, B. Prostaglandin Endoperoxides. Novel Transformations of Arachidonic Acid in Human Platelets. Proc. Natl. Acad. Sci. USA 1974, 71, 3400–3404. [CrossRef] [PubMed] 5. Ikei, K.N.; Yeung, J.; Apopa, P.L.; Ceja, J.; Vesci, J.; Holman, T.R.; Holinstat, M. Investigations of Human Platelet-Type 12- Lipoxygenase: Role of Lipoxygenase Products in Platelet Activation. J. Lipid Res. 2012, 53, 2546–2559. [CrossRef] [PubMed] 6. Ivanov, I.; Heydeck, D.; Hofheinz, K.; Roffeis, J.; O’Donnell, V.B.; Kuhn, H.; Walther, M. Molecular Enzymology of Lipoxygenases. Arch. Biochem. Biophys. 2010, 503, 161–174. [CrossRef] [PubMed] 7. Kulkarni, A.; Nadler, J.L.; Mirmira, R.G.; Casimiro, I. Regulation of Tissue Inflammation by 12-Lipoxygenases. Biomolecules 2021, 11, 717. [CrossRef] 8. Nugteren, D.H. Arachidonate Lipoxygenase in Blood Platelets. Biochim. Biophys. Acta (BBA)/Lipids Lipid Metab. 1975, 380, 299–307. [CrossRef] 9. Bürger, F.; Krieg, P.; Marks, F.; Fürstenberger, G. Positional- and Stereo-Selectivity of Fatty Acid Oxygenation Catalysed by Mouse (12S)-Lipoxygenase Isoenzymes. Biochem. J. 2000, 348, 329–335. [CrossRef] 10. Yu, J.Y.; Lee, J.J.; Jung, J.K.; Min, Y.K.; Ma, J.Y.; Kim, T.J.; Lee, M.Y.; Yun, Y.P. Anti-Platelet Activity of Diacetylated Obovatol through Regulating Cyclooxygenase and Lipoxygenase Activities. Arch. Pharm. Res. 2012, 35, 2191–2198. [CrossRef] 11. Katoh, A.; Ikeda, H.; Murohara, T.; Haramaki, N.; Ito, H.; Imaizumi, T. Platelet-Derived 12-Hydroxyeicosatetraenoic Acid Plays an Important Role in Mediating Canine Coronary Thrombosis by Regulating Platelet Glycoprotein IIb/IIIa Activation. Circulation 1998, 98, 2891–2898. [CrossRef] 12. Yoshimoto, T.; Takahashi, Y. Arachidonate 12-Lipoxygenases. Prostaglandins Other Lipid Mediat. 2002, 68–69, 245–262. [CrossRef] 13. Matsuda, S.; Murakami, J.; Yamamoto, Y.; Konishi, Y.; Yokoyama, C.; Yoshimoto, T.; Yamamoto, S.; Mimura, Y.; Okuma, M. Decreased Messenger RNA of Arachidonate 12-Lipoxygenase in Platelets of Patients with Myeloproliferative Disorders. BBA-Mol. Basis Dis. 1993, 1180, 243–249. [CrossRef] [PubMed] Int. J. Mol. Sci. 2023, 24, 6064 27 of 29 14. Sekiya, F.; Takagi, J.; Sasaki, K.; Kawajiri, K.; Kobayashi, Y.; Sato, F.; Saito, Y. Feedback Regulation of Platelet Function by 12 S-Hydroxyeicosatetraenoic Acid: Inhibition of Arachidonic Acid Liberation from Phospholipids. Biochim. Biophys. Acta (BBA)/Lipids Lipid Metab. 1990, 1044, 165–168. [CrossRef] 15. Holm, A.C.B.S.; Grenegård, M.; Öllinger, K.; Lindström, E.G. Inhibition of 12-Lipoxygenase Reduces Platelet Activation and Prevents Their Mitogenic Function. Platelets 2014, 25, 111–117. [CrossRef] 16. Adili, R.; Tourdot, B.E.; Mast, K.; Yeung, J.; Freedman, J.C.; Green, A.; Luci, D.K.; Jadhav, A.; Simeonov, A.; Maloney, D.J.; et al. First Selective 12-LOX Inhibitor, ML355, Impairs Thrombus Formation and Vessel Occlusion in Vivo with Minimal Effects on Hemostasis. Arterioscler. Thromb. Vasc. Biol. 2017, 37, 1828–1839. [CrossRef] [PubMed] 17. Tsai, W.C.; Aleem, A.M.; Tena, J.; Rivera-Velazquez, M.; Brah, H.S.; Tripathi, S.; D’silva, M.; Nadler, J.L.; Kalyanaraman, C.; Jacobson, M.P.; et al. Docking and Mutagenesis Studies Lead to Improved Inhibitor Development of ML355 for Human Platelet 12-Lipoxygenase. Bioorganic Med. Chem. 2021, 46, 116347. [CrossRef] 18. Contursi, A.; Tacconelli, S.; Hofling, U.; Bruno, A.; Dovizio, M.; Ballerini, P.; Patrignani, P. Biology and Pharmacology of Platelet-Type 12-Lipoxygenase in Platelets, Cancer Cells, and Their Crosstalk. Biochem. Pharmacol. 2022, 205, 115252. [CrossRef] 19. Dobrian, A.D.; Lieb, D.C.; Cole, B.K.; Taylor-Fishwick, D.A.; Chakrabarti, S.K.; Nadler, J.L. Functional and Pathological Roles of the 12- and 15-Lipoxygenases. Prog. Lipid Res. 2011, 50, 115–131. [CrossRef] 20. Serhan, C.N.; Romano, M. Lipoxin Biosynthesis and Actions: Role of the Human Platelet LX-Synthase. J. Lipid Mediat. Cell Signal. 1995, 12, 293–306. [CrossRef] 21. Kutzner, L.; Goloshchapova, K.; Heydeck, D.; Stehling, S.; Kuhn, H.; Schebb, N.H. Mammalian ALOX15 Orthologs Exhibit Pronounced Dual Positional Specificity with Docosahexaenoic Acid. Biochim. Biophys. Acta—Mol. Cell Biol. Lipids 2017, 1862, 666–675. [CrossRef] 22. Aleem, A.M.; Tsai, W.C.; Tena, J.; Alvarez, G.; Deschamps, J.; Kalyanaraman, C.; Jacobson, M.P.; Holman, T. Probing the Electrostatic and Steric Requirements for Substrate Binding in Human Platelet-Type 12-Lipoxygenase. Biochemistry 2019, 58, 848–857. [CrossRef] 23. Freedman, C.; Tran, A.; Tourdot, B.E.; Kalyanaraman, C.; Perry, S.; Holinstat, M.; Jacobson, M.P.; Holman, T.R. Biosynthesis of the Maresin Intermediate, 13S,14S-Epoxy-DHA, by Human 15-Lipoxygenase and 12-Lipoxygenase and Its Regulation through Negative Allosteric Modulators. Biochemistry 2020, 59, 1832–1844. [CrossRef] 24. Saura, P.; Suardíaz, R.; Masgrau, L.; Lluch, J.M.; González-Lafont, À. Unraveling How Enzymes Can Use Bulky Residues to Drive Site-Selective c-h Activation: The Case of Mammalian Lipoxygenases Catalyzing Arachidonic Acid Oxidation. ACS Catal. 2014, 4, 4351–4363. [CrossRef] 25. Weinberg, D.R.; Gagliardi, C.J.; Hull, J.F.; Murphy, C.F.; Kent, C.A.; Westlake, B.C.; Paul, A.; Ess, D.H.; McCafferty, D.G.; Meyer, T.J. Proton-Coupled Electron Transfer. Chem. Rev. 2012, 112, 4016–4093. [CrossRef] [PubMed] 26. Perry, S.C.; Kalyanaraman, C.; Tourdot, B.E.; Conrad, W.S.; Akinkugbe, O.; Freedman, J.C.; Holinstat, M.; Jacobson, M.P.; Holman, T.R. 15-Lipoxygenase-1 Biosynthesis of 7S,14S-DiHDHA Implicates 15-Lipoxygenase-2 in Biosynthesis of Resolvin D5. J. Lipid Res. 2020, 61, 1087–1103. [CrossRef] [PubMed] 27. Yamaguchi, A.; Stanger, L.; Freedman, C.J.; Standley, M.; Hoang, T.; Adili, R.; Tsai, W.C.; van Hoorebeke, C.; Holman, T.R.; Holinstat, M. DHA 12-LOX-Derived Oxylipins Regulate Platelet Activation and Thrombus Formation through a PKA-Dependent Signaling Pathway. J. Thromb. Haemost. 2021, 19, 839–851. [CrossRef] 28. Dalli, J.; Zhu, M.; Vlasenko, N.A.; Deng, B.; Haeggström, J.Z.; Petasis, N.A.; Serhan, C.N. The Novel 13S,14S-Epoxy-Maresinis Converted by Human Macrophages to Maresin 1 (MaR1), Inhibits Leukotriene A4 Hydrolase (LTA4H), and Shifts Macrophage Phenotype. FASEB J. 2013, 27, 2573–2583. [CrossRef] 29. Tang, S.; Wan, M.; Huang, W.; Stanton, R.C.; Xu, Y. Maresins: Specialized Proresolving Lipid Mediators and Their Potential Role in Inflammatory-Related Diseases. Mediat. Inflamm. 2018, 2018, 2380319. [CrossRef] 30. Serhan, C.N.; Dalli, J.; Karamnov, S.; Choi, A.; Park, C.-K.; Xu, Z.-Z.; Ji, R.-R.; Zhu, M.; Petasis, N.A. Macrophage Proresolving Mediator Maresin 1 Stimulates Tissue Regeneration and Controls Pain. FASEB J. 2012, 26, 1755–1765. [CrossRef] 31. Vogel, R.; Jansen, C.; Roffeis, J.; Reddanna, P.; Forsell, P.; Claesson, H.E.; Kuhn, H.; Walther, M. Applicability of the Triad Concept for the Positional Specificity of Mammalian Lipoxygenases. J. Biol. Chem. 2010, 285, 5369–5376. [CrossRef] 32. Chen, X.-S.; Funk, C.D. Structure-function Properties of Human Platelet 12-lipoxygenase: Chimeric Enzyme and in Vitro Mutagenesis Studies. FASEB J. 1993, 7, 694–701. [CrossRef] [PubMed] 33. Neau, D.B.; Gilbert, N.C.; Bartlett, S.G.; Boeglin, W.; Brash, A.R.; Newcomer, M.E. The 1.85 Å Structure of an 8R-Lipoxygenase Suggests a General Model for Lipoxygenase Product Specificity. Biochemistry 2009, 48, 7906–7915. [CrossRef] 34. Neau, D.B.; Bender, G.; Boeglin, W.E.; Bartlett, S.G.; Brash, A.R.; Newcomer, M.E. Crystal Structure of a Lipoxygenase in Complex with Substrate: The Arachidonic Acid-Binding Site of 8R-Lipoxygenase. J. Biol. Chem. 2014, 289, 31905–31913. [CrossRef] [PubMed] 35. Newcomer, M.E.; Brash, A.R. The Structural Basis for Specificity in Lipoxygenase Catalysis. Protein Sci. 2015, 24, 298–309. [CrossRef] [PubMed] 36. Xu, S.; Mueser, T.C.; Marnett, L.J.; Funk, M.O. Crystal Structure of 12-Lipoxygenase Catalytic-Domain-Inhibitor Complex Identifies a Substrate-Binding Channel for Catalysis. Structure 2012, 20, 1490–1497. [CrossRef] [PubMed] 37. Yokoyama, C.; Shinjo, F.; Yoshimoto, T.; Yamamoto, S.; Oates, J.A.; Brash, A.R. Arachidonate 12-Lipoxygenase Purified from Porcine Leukocytes by Immunoaffinity Chromatography and Its Reactivity with Hydroperoxyeicosatetraenoic Acids. J. Biol. Chem. 1986, 261, 16714–16721. [CrossRef] [PubMed] Int. J. Mol. Sci. 2023, 24, 6064 28 of 29 38. Ueda, N.; Hiroshima, A.; Natsui, K.; Shinjo, F.; Yoshimoto, T.; Yamamoto, S.; Ii, K.; Gerozissis, K.; Dray, F. Localization of Arachidonate 12-Lipoxygenase in Parenchymal Cells of Porcine Anterior Pituitary. J. Biol. Chem. 1990, 265, 2311–2316. [CrossRef] 39. Reddy, R.G.; Yoshimoto, T.; Yamamoto, S.; Marnett, L.J. Expression, Purification, and Characterization of Porcine Leukocyte 12-Lipoxygenase Produced in the Methylotrophic Yeast, Pichia Pastoris. Biochem. Biophys. Res. Commun. 1994, 205, 381–388. [CrossRef] 40. Suzuki, H.; Kishimoto, K.; Yoshimoto, T.; Yamamoto, S.; Kanai, F.; Ebina, Y.; Miyatake, A.; Tanabe, T. Site-Directed Mutagenesis Studies on the Iron-Binding Domain and the Determinant for the Substrate Oxygenation Site of Porcine Leukocyte Arachidonate 12-Lipoxygenase. Biochim. Biophys. Acta (BBA)/Lipids Lipid Metab. 1994, 1210, 308–316. [CrossRef] 41. Saura, P.; Kaganer, I.; Heydeck, D.; Lluch, J.M.; Kühn, H.; González-Lafont, À. Mutagenesis of Sequence Determinants of Truncated Porcine ALOX15 Induces Changes in the Reaction Specificity by Altering the Catalytic Mechanism of Initial Hydrogen Abstraction. Chem.—Eur. J. 2018, 24, 962–973. [CrossRef] 42. Jumper, J.; Evans, R.; Pritzel, A.; Green, T.; Figurnov, M.; Ronneberger, O.; Tunyasuvunakool, K.; Bates, R.; Žídek, A.; Potapenko, A.; et al. Highly Accurate Protein Structure Prediction with AlphaFold. Nature 2021, 596, 583–589. [CrossRef] [PubMed] 43. Varadi, M.; Anyango, S.; Deshpande, M.; Nair, S.; Natassia, C.; Yordanova, G.; Yuan, D.; Stroe, O.; Wood, G.; Laydon, A.; et al. AlphaFold Protein Structure Database: Massively Expanding the Structural Coverage of Protein-Sequence Space with High-Accuracy Models. Nucleic Acids Res. 2022, 50, D439–D444. [CrossRef] 44. Aleem, A.M.; Jankun, J.; Dignam, J.D.; Walther, M.; Kühn, H.; Svergun, D.I.; Skrzypczak-Jankun, E. Human Platelet 12- Lipoxygenase, New Findings about Its Activity, Membrane Binding and Low-Resolution Structure. J. Mol. Biol. 2008, 376, 193–209. [CrossRef] [PubMed] 45. Shang, W.; Ivanov, I.; Svergun, D.I.; Borbulevych, O.Y.; Aleem, A.M.; Stehling, S.; Jankun, J.; Kühn, H.; Skrzypczak-Jankun, E. Probing Dimerization and Structural Flexibility of Mammalian Lipoxygenases by Small-Angle X-Ray Scattering. J. Mol. Biol. 2011, 409, 654–668. [CrossRef] [PubMed] 46. Romero-Téllez, S.; Cruz, A.; Masgrau, L.; González-Lafont, À.; Lluch, J.M. Accounting for the Instantaneous Disorder in the Enzyme-Substrate Michaelis Complex to Calculate the Gibbs Free Energy Barrier of an Enzyme Reaction. Phys. Chem. Chem. Phys. 2021, 23, 13042–13054. [CrossRef] 47. Martínez-Rosell, G.; Giorgino, T.; De Fabritiis, G. PlayMolecule ProteinPrepare: A Web Application for Protein Preparation for Molecular Dynamics Simulations. J. Chem. Inf. Model. 2017, 57, 1511–1516. [CrossRef] 48. Jones, G.; Willett, P.; Glen, R.C.; Leach, A.R.; Taylor, R. Development and Validation of a Genetic Algorithm for Flexible Docking. J. Mol. Biol. 1997, 267, 727–748. [CrossRef] 49. Tian, C.; Kasavajhala, K.; Belfon, K.A.A.; Raguette, L.; Huang, H.; Migues, A.N.; Bickel, J.; Wang, Y.; Pincay, J.; Wu, Q.; et al. Ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution. J. Chem. Theory Comput. 2020, 16, 528–552. [CrossRef] 50. Izadi, S.; Anandakrishnan, R.; Onufriev, A.V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5, 3863–3871. [CrossRef] 51. 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] 52. 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] 53. 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] 54. 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. Gaussian 16; Gaussian, Inc.: Wallingford, CT, USA, 2016. 55. Li, P.; Merz, K.M. MCPB.Py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016, 56, 599–604. [CrossRef] 56. Seminario, J.M. Calculation of Intramolecular Force Fields from Second-Derivative Tensors. Int. J. Quantum Chem. 1996, 60, 1271–1277. [CrossRef] 57. 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 20; University of California: San Francisco, CA, USA, 2020. 58. Case, D.A.; Betz, R.M.; Cerutti, D.S.; Cheatham III, T.E.; Darden, T.A.; Duke, R.E.; Giese, T.J.; Gohlke, H.; Goetz, A.W.; Homeyer, N.; et al. Amber 16; University of California: San Francisco, CA, USA, 2016. 59. 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] 60. Salomon-Ferrer, R.; Götz, A.W.; Poole, D.; Le Grand, S.; Walker, R.C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878–3888. [CrossRef] [PubMed] 61. Leach, A.R. Molecular Modelling. Principles and Applications, 2nd ed.; Pearson Education: London, UK, 2001; ISBN 978-0-582-38210-7. 62. 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] 63. 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] Int. J. Mol. Sci. 2023, 24, 6064 29 of 29 64. Canyelles-Niño, M. MolBioMedUAB/RCBS.Py: 0.2.0bis. GitHub: San Franciso, CA, USA, 2023. [CrossRef] 65. Hunter, J.D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007, 9, 90–95. [CrossRef] 66. 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] 67. Metz, S.; Kästner, J.; Sokol, A.A.; Keal, T.W.; Sherwood, P. ChemShell-a Modular Software Package for QM/MM Simulations. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 101–110. [CrossRef] 68. 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] 69. 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] 70. Becke, A.D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652. [CrossRef] 71. Hariharan, P.C.; Pople, J.A. The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213–222. [CrossRef] 72. Hay, P.J.; Wadt, W.R. Ab Initio Effective Core Potentials for Molecular Calculations. Potentials for the Transition Metal Atoms Sc to Hg. J. Chem. Phys. 1985, 82, 270–283. [CrossRef] 73. 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] 74. Liu, D.C.; Nocedal, J. On the Limited Memory BFGS Method for Large Scale Optimization. Math. Program. 1989, 45, 503–528. [CrossRef] 75. 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] 76. 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] 77. Heyden, A.; Bell, A.T.; Keil, F.J. Efficient Methods for Finding Transition States in Chemical Reactions: Comparison of Improved Dimer Method and Partitioned Rational Function Optimization Method. J. Chem. Phys. 2005, 123, 224101. [CrossRef] [PubMed] 78. 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] 79. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera—A Visualiza- tion System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605–1612. [CrossRef] [PubMed] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.