This is the accepted version of the journal article: Cebrián Prats, Anna; Pinto, Alexandre; González Lafont, Àngels; [et al.]. «The role of acetylated cyclooxygenase-2 in the biosynthesis of resolvin precursors derived from eicosapentaenoic acid». Organic & Biomolecular Chemistry, Vol. 20, Issue 6 (February 2022), p. 1260-1274. DOI 10.1039/D1OB01932E This version is available at https://ddd.uab.cat/record/259977 under the terms of the license The Role of Acetylated Cyclooxygenase-2 in the Biosynthesis of Resolvin Precursors Derived from Eicosapentaenoic Acid Anna Cebrián-Prats1, Alexandre Pinto 2, Àngels González-Lafont*1,3, Pedro A. Fernandes*2, and José M. Lluch1,3 1Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain 2LAQV-Requimte, Faculty of Sciences, University of Porto, Rua do Campo Alegre S/N, Porto, Portugal. 3Institut de Biotecnologia i de Biomedicina (IBB), Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain Abstract Specialized pro-resolving lipid mediators (SPMs) are natural bioactive agents actively involved in inflammation resolution. SPMs act when uncontrolled inflammatory processes are developed, for instance, in patients of COVID-19 or other diseases. The so-called resolution pharmacology aims at developing new treatments based on the use of SPMs as agonists, which promote inflammation resolution without unwanted side effects. It has been shown that the biosynthesis of the SPMs called eicosapentaenoic acid (EPA)-derived E-series resolvins is initiated by aspirin-acetylated COX-2 from EPA, leading to 18-hydroperoxy-eicosapentaenoic acid (18-HpEPE). However, there are many open questions concerning the intriguing role of aspirin in the molecular mechanism of resolvins formation. Our MD simulations, combined with QM/MM calculations, show that the potential energy barriers for the H16-abstraction from EPA, required for forming 18-HpEPE, are higher than for the H13-abstraction, so explaining why 18-HpEPE is a marginal product of COX-2 catalysis. By contrast, in the aspirin- acetylated COX-2/EPA complex, the H16proS-abstraction energy barriers are somewhat lower than the H13proS energy barriers and much smaller than the H16-transfer barriers in the wild type COX-2/EPA system. Those results agree with the experimental observation that aspirin favours the synthesis of several SPMs known as aspirin-triggered resolvins. In the following step of the catalytic mechanism, the calculated O2 addition to C18 is preferred versus the addition to C14 what also agrees with 18R- HEPE and 18S-HEPE being the main products from EPA in aspirin-acetylated COX-2. 1 1. Introduction The acute and uncontrolled inflammation that manifests in many different diseases such as cancer, neurological, cardiovascular, and respiratory disorders, is also a clear pathology in Coronavirus Disease 2019 (COVID-19).1,2 COVID-19 patients are suffering from the hyperactivated inflammatory response, triggered by SARS-CoV-2, of their autoimmune system. This generalized inflammation sometimes produces multi-organ dysfunctions and death. For that reason, the study of inflammation resolution is now of public health concern. The mechanisms involved in those systemic inflammatory effects and the way to treat them has become the target of many undergoing biomedical studies.3 The resolution of inflammation is now recognized as an active process carried out by a superfamily of bioactive agents called specialized proresolvin lipid mediators (SPMs).4,5 SPMs are synthesized in macrophages from polyunsaturated fatty acids (PUFAs). An active process means here that the inflammatory agents (like leukotrienes or cytokines) do not simply dissipate when resolution takes place. The presence of SPMs is needed to stimulate inflammation resolution. SPMs inhibit those proinflammatory molecules and finally regenerate tissues to their normal structure by homeostasis. In humans, SPMs include three main types of mediators derived from ω-3 PUFAs. Those SPMs are named protectins, resolvins and, more recently, a new class of SPMs called maresins have been discovered.6,7 Protectins, D-resolvins and maresins are bioactive SPMs that come from docosahexaenoic acid (DHA). E-resolvins are generated from eicosapentaenoic acid (EPA).8 Both acids are long-chain essential PUFAs present in fish-oil and widely used as dietary supplements.9,10 DHA has 22 carbon atoms and six unsaturations, and EPA has 20 carbon atoms and five unsaturations. EPA is a substrate of cyclooxygenase-2 (COX-2) that catalyzes its conversion to prostaglandin PGH3 following, in principle, the same all-radical mechanism than for producing prostaglandin PGH2 from arachidonic acid (AA).5 However, EPA is a poorer substrate of COX-2 than AA and it is oxygenated at about 30% of the rate of AA.11 Also, the product profiles are different when comparing both substrates. The proportion of PGH3 versus hydroxy-eicosapentaenoic acids (HEPEs) secondary 2 products is around 50% whereas for PGH2 this proportion reaches 90%.12 The main secondary product of EPA oxidation that has been reported is 11R-hydroxy-eicosapentaenoic acid (11R-HEPE), followed by 15R-hydroxy-eicosapentaenoic acid (15R-HEPE) and a smaller amount of 18R-hydroxy- eicosapentaenoic acid (18R-HEPE).13 In some other experiments 14-hydroxy-eicosapentaenoic acid (14-HEPE) was included as a secondary product instead of 18-HEPE, but the characterization of 14- HEPE was not confirmed by using an authentic standard.12 Aspirin has a positive effect on the synthesis of several SPMs, whereas other drugs abolish timely resolution. In the case of EPA-derived E-series resolvins, their biosynthesis from EPA, catalyzed by aspirin-acetylated COX-2 (aspirin-triggered resolvins) (see Figure 1), begins by the production of 18- hydroperoxy-eicosapentaenoic acid (18-HpEPE) that it is subsequently reduced by a peroxidase to 18- HEPE. In a second lipoxygenation catalyzed by 5-LOX, a hydroperoxide intermediate is formed that evolves to an epoxide intermediate in the same 5-LOX active site, and that finally leads to Resolvin E1 (RvE1) via the action of a hydrolase. Resolvin E2 (RvE2) is formed by the reduction of the hydroperoxide intermediate via the action of a peroxidase. The biosynthetic pathway leading to the formation of E-resolvins, as well as the contribution of aspirin and the different enzymes involved in the mechanism, was elucidated by Serhan and coworkers.8,14 Both, RvE1 and RvE2, were initially isolated from in vitro and in vivo models of resolution,13,15 and the authors paid special attention to properly assign the regiospecificity and stereospecificity of the different products. Because resolvins were identified in mouse exudates, it was essential to establish their biosynthesis by human leukocytes and in human tissues.13,15 Therefore, using mass spectrometry, RvE1 and RvE2 were identified in the peripheral blood of healthy volunteers, some of whom had been given EPA supplements.8 The potent actions of RvE1 and RvE2 following their complete stereochemical assignments have been confirmed by many independent laboratories.7 3 Figure 1. Biosynthetic pathways to obtain Resolvin E1/E2, and 18S-Resolvin E1/E2 proposed by Serhan and coworkers.16,17 In the original experiments, only the 18R stereochemistry was assigned to the HEPE precursors of RvE1 and RvE2. However, in later experimental analysis, the chiral LC-MS/MS based lipidomics technique was applied and the 18S-HEPE isomer was also identified in healthy humans after the intake of aspirin and EPA.16,18 Serhan et al. explain that according to their data, aspirin increases the production of both 18R-and 18S-HEPE, the first one being still the major product. Then they conclude that the formation of 18S- and 18R-HEPE is catalyzed by aspirin-acetylated COX-2. The resolvins generated from 18S-HEPE are denominated as 18S-Resolvin E1 and 18S-Resolvin E2 as shown in Figure 1. It has been stated that the traditional therapy of inflammatory effects based on the intake of non- steroidal anti-inflammatory drugs (NSAIDs) that act as cyclooxygenase (COX) inhibitors must be 4 modified.18 Those drugs function by blocking the spread of inflammation and infection in the body. However, all of them are known to present unwanted side-effects and immune suppression. Aspirin jump-starts resolution because it triggers the biosynthesis of SPMs (this provides additional basis for the health benefits of aspirin). In this respect, Serhan has recently affirmed that the new therapy of the 21st century must be focused on promoting natural and programmed resolution processes by using SPMs as agonists.18 This new approach referred as resolution pharmacology is based on the idea that resolution of inflammation needs a certain degree of inflammation. That is, some doses of inflammation are capable to trigger the initiation of the resolution process so maintaining the delicate balance between inflammatory and anti-inflammatory agents without causing immunosuppression. In this study, we employed molecular dynamics simulations and QM/MM calculations to investigate the molecular details underpinning the catalytic mechanism of aspirin-acetylated COX-2 for the formation of EPA-derived products that are precursors of aspirin-triggered SPMs. Some differences with the same catalytic mechanism in COX-2 are discussed using a QM/MM analysis of the first chemical step (H-abstraction) of the two catalytic processes. Also, structural differences between the EPA binding modes within the active site of COX-2 and within the pocket of aspirin-acetylated COX-2 are highlighted. The results confirm the energetic preference in aspirin-acetylated COX-2 for the biosynthetic pathway leading to 18-HEPE versus the pathways leading to the secondary products 11R- HEPE, 15R-HEPE, and 14R-HEPE. 2. Methods 2.1 Model setup and parameterization In the model for the COX-2/EPA complex, the enzyme was built based on the murine COX-2/EPA coordinates (PDB code 3HS6).19 In this structure, the protoporphyrin IX group is coordinated to a Co3+ ion. The homodimer presents one EPA ligand bound to each monomer in a different conformation. We replaced the Co3+ ion by a Fe3+ one, and the heme prosthetic group was tethered to the axial 5 histidine residue (His388) to prevent the dissociation of the complex during the simulation. As in the COX-2/AA complex,19 only the monomer B exhibits the productive configuration for EPA in the cyclooxygenase channel, so this binding pose was used for this work. The AMBER ff14SB force field20 was employed to define the protein residues. The missing hydrogen atoms were added using the leap module of AMBER16.21 All the residues are considered in their standard protonation states assuming a pH = 7. The heme group and the Fe3+ ion, residue His388, and the tyrosyl radical were parameterized, as explained in our previous studies. That is, the parameters of the heme group were adopted from a recent study22 by another group even though we recalculated the partial charges of the penta- coordinated high spin (Fe3+) group employing the RESP algorithm23 at the HF/6-31G(d,p) level for the QM calculations. The GAFF24 parameters were used for His388. For the tyrosyl residue, we developed specific force constants for bonds and angles, fitting the stretching and bending energy profiles of this residue at the B3LYP/6-31G(d,p) level. For this residue, the RESP charges were also derived at the B3LYP/6-31G(d) level using the optimized geometry at the same level. As for the van der Waals parameters, we used GAFF values. Finally, the parameters for the EPA substrate and the C11-C15 and C15-C18 pentadienyl radicals were taken from the GAFF forcefield of AMBER16. The total system was solvated in a nearly cubic box (99.8 Å x 110.9 Å x 97.2 Å) with TIP3P25 water molecules. A distance of 15 Å was established between the atoms of the enzyme and the edge of the box. Water molecules closer than 2.2 Å to any atom of the EPA ligand or the enzyme were removed. The total charge of the system was neutralized by including five Na+ ions. The final solvated COX-2/EPA Michaelis complex has 81558 atoms. The model for the aspirin-acetylated COX-2/EPA Michaelis complex was built based on monomer B of the aspirin-acetylated COX-2/AA complex used in our previous study and replacing AA by EPA.26 The same force field described in the previous paragraph was used for the protein, His388, the tyrosyl radical, and the EPA substrate. The parameters for the aspirin-acetylated Ser530 were taken from our previous parameterization26 for the aspirin-acetylated COX-2/AA system that was based on the crystal structure of the human aspirin-acetylated COX-2 (PDB code 5F19).27 The total system was also solvated 6 with TIP3P water molecules in a box with the same dimensions than for the COX-2/EPA system. The final solvated aspirin-acetylated COX-2/EPA Michaelis complex has 81665 atoms. 2.2 MD Simulations The resulting COX-2/EPA complex was optimized in three steps. First, we applied harmonic restraints on the substrate, the enzyme and the heme prosthetic group, keeping the solvation waters free. Secondly, only the substrate coordinates and the heme prosthetic group have been restrained. Finally, we minimized the whole system, only restraining the protein side chains. Afterwards, we have carried out the MD simulations under periodic boundary conditions using the same protocol as in our study of the COX-2/AA system.28 The interatomic interactions are evaluated with a cutoff of 10 Å for all Lennard-Jones and electrostatic contributions and using the particle-mesh Ewald method29 to treat long-range electrostatic effects. A heating step of 200 ps using Langevin dynamics30 with weak restraints on the substrate and the heme group was carried out to take the complex from 0 to 300 K in at constant volume. Then, an NPT MD trajectory of 200 ps (using an isotropic weak-coupling algorithm and the Berendsen barostat31 with weak restraints on the protein side chains) was performed until the density of the system converged around a value of 1 g/cm3. Finally, we performed an MD simulation of 10 ns to equilibrate the Michaelis complex and a production MD trajectory of 100 ns without restraints at 300 K within an NVT ensemble. A time step of 2 fs has been chosen. We repeated this protocol for two MD simulations: one with the Tyr385 residue and another with the Tyr385 radical. The simulations were performed using the AMBER 16 GPU (CUDA) version of the PMEMD package.32,33 The resulting model for the aspirin-acetylated COX-2/EPA system was minimized according to the three-stages optimization described above. In this case, we carried out one MD simulation with the tyrosyl radical under periodic boundary conditions for the aspirin-acetylated COX-2/EPA complex, applying the same protocol than for the COX-2/EPA complex. 7 Finally, representative snapshots (see below) of the MD simulations were used to perform the QM/MM calculations. For the COX-2/EPA system, those snapshots were selected to initiate the H- abstraction reaction which is the first step of the all-radical mechanism. For the aspirin-acetylated COX-2/EPA complex the selected snapshots were ready to initiate the H-abstraction process and then were also used to study the addition of O2, the second step in the reaction pathway of 18-HEPE formation. 2.3 QM/MM calculations The QM/MM calculations were carried out using the ONIOM (QM:MM) approach34 as implemented in the Gaussian09 software.35 The B3LYP density functional36,37,38 with the 6-31G(d,p) basis set was used in the QM layer and the amber force field39 was used in the MM layer as implemented in Gaussian09. The interaction between the QM and MM layers was treated with the electrostatic embedding scheme. Hydrogen link atoms were used to fill the vacant valences caused by truncating bonds across the QM and MM layers. The QM layer included the side chain of Tyr385 radical up to the beta-carbon and all atoms of the EPA substrate from C10 to C20 for the H-abstraction processes (a total of 42 atoms). In the QM/MM calculations for the oxygen addition processes, the oxygen molecule is added to the QM region (a total of 44 atoms) (see Figure 2). A shell of solvent was kept around the protein, which included all water molecules within 3 Å of any protein atom. During geometry optimizations, all MM atoms within a radius of 15 Å of the QM layer were kept free and the remaining atoms were frozen. This process avoids drifting along multiple minima without affecting the accuracy of the results.40 The total charge of the QM layer was zero and the spin multiplicity was 2. 8 Figure 2. QM region depicted in blue for the QM/MM calculations of the H-abstraction processes. Most H atoms in the QM region have not been shown for the sake of clarity. The red wavy lines indicate the location of the two link atoms needed at the boundary between the QM and MM regions. For the oxygen addition processes, the oxygen molecule has been added to the QM region. Linear transit scans along trial reaction coordinates were performed to generate trial transition structures. The latter were subsequently fully optimized, and their nature was confirmed by means of the analysis of the number of imaginary frequencies (only one for transition state structures). The corresponding reactants and products were also fully optimized. For the oxygen additions the energy minima of the energy profiles were fully optimized whereas the maxima of the energy profiles defined as the corresponding approximate transition state structures, were used without further optimization. QM/MM single point energy calculations in the stationary states was made at the more accurate B3LYP/6-311+G(2d,2p):AMBER level of theory; The D3 version of Grimme's dispersion with Becke-Johnson damping dispersion corrections41 was added afterwards. 3. Results and Discussion 3.1 Molecular dynamics simulations of the COX-2/EPA complex We first carried out the simulation protocol to equilibrate the wild-type COX-2/EPA complex. Two different systems were modeled: one with EPA in the active site in the presence of Tyr385, and a second one with EPA and the Tyr385-O· radical. As already indicated, the peroxidase cycle takes place 9 independently of the cyclooxygenase activity,8,14 so it is not unfeasible that EPA could enter the cyclooxygenase active site before or after the tyrosyl radical of Tyr385 is formed. The behavior of the root mean-square deviation (RMSD) of the protein α-carbons with respect to the first structure shows that the protein is equilibrated in both systems, in particular if the N-terminal region is not considered (see Figure S1). Also, the most flexible part of the protein is its N-terminal region which shows a rather big displacement and fluctuations (see the RMSF plot in Figure S2). As indicated above, the first step of the overall mechanism consists of a H-transfer from the EPA substrate to the Tyr385-O· radical. Considering that the COX-2/EPA produces PGH3 as a major product, and 11-HEPE, 15-HEPE 14-HEPE and 18-HEPE as secondary products, the hydrogen abstraction could be at C13 as well as C16. To establish the viability of this first chemical reaction, we analyzed the evolution of several relevant interatomic distances along the two MD trajectories. Thus, we recorded the H13proX-OTyr385 and H16proX-OTyr385 distances (where proX stands for either proS or proR hydrogen) along both trajectories (see Figure 3). The H13proX abstraction will lead to PGH3, the main product of EPA oxidation in COX-2, and to two of its most abundant secondary products, 11- hydroperoxy-eicosapentaenoic acid and 15-hydroperoxy-eicosapentaenoic acid, whereas the H16proX abstraction will lead to 18-hydroperoxy-eicosapentaenoic acid, a minor product of EPA oxidation in COX-2 according to some previous experimental results.12,13 The comparison between the pre- catalytic structures for both H-abstraction processes will provide molecular clues on how the EPA binding mode at the active site might influence the regioselectivity shown by COX-2. We will define the pre-catalytic structures as those where at least one well-oriented hydrogen atom at C13 or C16 is closer than 3.0 Å from the oxygen atom of the FeIII−OH- cofactor, then these structures being ready to react. It can be observed in Figure 3 that the four distances H13proX-OTyr385 fluctuate along the two trajectories. During half simulation with Tyr385, H13proR remains closer to the tyrosine residue than H13proS, although both hydrogens interchange quite often (the average H13proR-OTyr385 and H13proS- OTyr385 distances are 3.27 Å and 4.08 Å, respectively). After a rearrangement that takes place after 55 ns, both H13 distances increase, and also the difference between them. In contrast, in the simulation 10 with the Tyr385-O· radical, the H13proR atom is closer to the tyrosine residue up to 55 ns (with an average H13proR-OTyr385 distance of 3.01 Å versus 4.37 Å for the H13proS-OTyr385 distance), but then the substrate also rearranges, and H13proR goes further away. In the last 45 ns, the H13proR and H13proS average distances are very similar, although H13proS presents more fluctuations that approach this atom to the Tyr385 radical. As for the H16proX-OTyr385 distances, it can be seen in Figure 3 that they show a somewhat reverse behaviour in comparison with the H13proX-OTyr385 distances. That is, in the simulation with Tyr385, the H16proR atom is closer to the tyrosine residue up to 55 ns, but then the substrate rearranges and H16proR goes further away. In the last 45 ns, the H16proR and H16proS average distances are very similar and both hydrogens remain too far away from the tyrosine residue to be abstracted. When the Tyr385-O· radical is present, H16proR remains in average closer to the tyrosine residue than H16proS (with an average H16proR-OTyr385 distance of 3.58 Å during the first 55 ns, versus 5.14 Å for the H16proS-OTyr385 distance), although both hydrogens interchange. During the last 45 ns both H16proR and H16proS remain in average at a similar distance from the tyrosyl radical. All the calculated average distances for the two MD simulations are collected in Table S1. 11 Figure 3. Evolution of the H13proX-OTyr385 and H16proX-OTyr385 distances along the two 100ns MD simulations of the COX-2/EPA Michaelis complex. On the right side, these distances for the MD simulation with the Tyr385-O· radical are represented. On the left side, the same distances for the MD simulation with Tyr385 are shown. To determine how many pre-catalytic structures appear along the simulations, one snapshot each 10 ps along both 100 ns MD simulations was saved. Then, a filtering of those structures was done according to the following geometrical conditions: d(HZproX-OTyr385) < 3 Å and d(HZproX-OTyr385)< d(CZ- OTyr385), to select the most adequate structures to initiate the all-radical mechanism. In Table 1 the percentage of pre-catalytic structures found for each simulation are given. It can be observed that the number of structures that accomplish the two geometrical conditions for H13proR and H16proR is bigger than for H13proS and H16proS, respectively, both for the simulation with Tyr385 and the simulation with the Tyr385 radical. Also, the structures that could initiate H13proR abstraction are more numerous than those corresponding to H16proR abstraction, especially in the trajectory with Tyr385. This result is 12 in agreement with PGH3 being the principal product12 oxidation by COX-2 as for that the all-radical mechanism has to be initiated by H13-abstraction. Our results also point to a more favorable HproR than HproS abstraction. Table 1. The percentage of feasible reactive structures (pre-catalytic structures) to initiate the all- radical mechanism by the abstraction of a hydrogen at C13 or C16 of EPA. The selection is performed along the two 100 ns MD simulations with either Tyr385 or the Tyr385-O· radical in the COX-2/EPA complex. H13proR H13proS H16proR H16proS Tyr385-OH 25.8 1.92 1.63 0 Tyr385-O· 32.5 1.23 16.8 6.9 We have seen in Table 1 that the number of pre-catalytic structures is clearly bigger in the case of the tyrosyl radical than with Tyr385. Then only the reactive structures of the simulation with the tyrosyl radical have been considered to initiate the calculation of the QM/MM H-abstraction energy profiles. On one hand, from 3252 pre-catalytic structures with the H13proR near the oxygen acceptor, four structures (snapshots I to IV) have been selected randomly. On the other hand, from 1677 pre-catalytic structures with the H16proR near the oxygen acceptor, three structures (snapshots V to VII) have been selected randomly. The initial HZproR-OTyr385 and CZ-OTyr385 distances for the seven chosen structures are given in Table 2. 13 Table 2. Initial distances (in Å) corresponding to the three atoms directly involved in the hydrogen abstraction (first step) of the all-radical mechanism for the seven selected snapshots in the case of the COX-2/EPA complex in the presence of the Tyr385-O· radical. H13proR H16proR Snapshots d(H-O) d(C-O) Snapshots d(H-O) d(C-O) I 2.51 3.53 V 2.34 3.34 II 2.45 3.50 II 2.85 3.66 III 2.95 3.87 VI 2.36 2.88 IV 2.65 3.56 VII 2.45 3.88 Snapshots I, II, III and IV were selected as adequate for the H13proR abstraction process, whereas snapshots V, VI and VII were considered appropiate for the H16proR transfer reaction. A similar approach to generate pre-catalytic structures was used for the COX-1/AA system.42 Snapshot II presents both hydrogens close enough to the Tyr385-O· radical to be a potential candidate for both hydrogen abstractions. The geometry of snapshots II, IV, V, and VI are compared in Figure 4. In the four snapshots, EPA presents a more “bulged L-shaped” binding mode than in the ones obtained in the simulations of the COX-2/AA complex28 due, in part, to a different conformation of the Ser530 side chain. At the opening of the channel, the EPA carboxylate interacts with the side chains of Arg513 and Arg120 in snapshots II and IV, and with Tyr355 and Arg120 in structures V and VI (see also Figure S3 and Table S2). The EPA ω-end extends along the hydrophobic groove above Ser530. The main overall difference between those snapshots is that the EPA ligand is less extended in those structures for which only the abstraction of H16proR could be feasible according to the initial H16proR-OTyr385 distance. 14 Figure 4. Binding modes of the EPA substrate inside the hydrophobic groove of COX-2 in the presence of the Tyr385-O· radical for snapshots II, IV, V and VI. 3.2 QM/MM calculations of the H13proR/H16proR abstraction processes in the COX-2/EPA complex As the first step in building up the H-abstraction energy profiles, we optimized the minimum energy structure for each reactant complex on the QM/MM potential energy surface at the B3LYP/6- 31G(d,p):AMBER level. Then, we calculated the potential energy profiles at discrete values of a reaction coordinate, defined as the distance corresponding to the forming bond (HzproR − OTyr385) for both the H13proR and H16proR atoms, from where we took the minima and maximum for further full optimization. Thus, we fully optimized the reactants, product and the transition state structures of the H-abstraction process. In the case of the two H16-abstractions, the energy barriers turn out to be quite high (see below), and the transition state structures were just approximated as the energy maxima of the corresponding potential energy profiles without further optimization. In Table 3, the B3LYP/6- 31G(d,p):AMBER potential energy barriers and the reaction energies are given for those snapshots 15 that have been able to reach the products (snapshots I, II, IV, and V). The calculation of the energy profiles starting from the minimum energy structures of snapshots III, VI, and VII did not succeed in leading to the hydrogen abstraction product. In the case of snapshot III the α-helix of the Ser530 residue is deformated, whereas in snapshot VI the deformation of the transition state structure geometry disables the calculation of the remaining potential energy points to reach the product. In turn, optimization of snapshot VII moves both the H13proR and H16proR hydrogen atoms far from the oxygen acceptor of the tyrosyl radical. The B3LYP/6-31G(d,p):AMBER potential energy barriers for H13proR abstraction are clearly lower than for the H16proR abstraction. These results confirm that the H13-transfer reaction is more favourable versus the H16-abstraction from EPA in COX-2, in agreement with the all-radical mechanism that initiates by a H13-abstraction step. It is also worth pointing out that the H13-abstraction energy barriers for EPA are, on average, higher than those for AA,28 also in agreement with the experimental results that have proved that EPA is a poorer substrate of COX-2 than AA.11 The lowest calculated barrier for the H13 abstraction of AA was 16.4 kcal/mol, whereas the lowest barrier for EPA is 22.3 kcal/mol. When we analyzed the QM and MM energy contributions to those different potential energy barriers (see Figure S4), we obtained that the QM contribution is by far the determinant factor. Along the reaction coordinates for snapshots I, II, IV, and V the MM energy contribution is small, and it is nearly constant as the H-transfer advances, showing only a slightly change after the transition state region. Table 3. Potential energy barriers and reaction energies (in kcal/mol) at the B3LYP/6-31G(d,p):AMBER level corresponding to the HZproR hydrogen abstractions for those snapshots that have been able to reach the products in the COX-2/EPA complex, as well as the corresponding values obtained by means of QM/MM single point energy calculations at the B3LYP/6-311+G(2d,2p):AMBER//B3LYP/6- 31G(d,p):AMBER, and including a posteriori GD3BJ dispersion corrections. 16 B3LYP/6-311+G B3LYP/6-311+G Snapshots B3LYP/6-31G(d,p) (2d,2p) (2d,2p) - GD3BJ H13proR ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E I 25.0 -1.0 27.5 -1.5 25.1 -2.8 II 26.9 -2.5 27.7 -3.5 25.0 -5.4 IV 22.3 -6.9 22.8 -8.7 19.2 -11.5 H16proR ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E II 31.9 2.5 33.5 2.0 34.6 2.5 V 29.1 -1.9 29.9 -3.2 27.9 -5.5 As seen in Table 3 the QM/MM single point energy results at the B3LYP/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER level qualitatively agree with the fully optimized B3LYP/6-31G(d,p):AMBER ones. The inclusion of the dispersion corrections lowers most of the potential energy barriers,43,44 but the relative order of those barriers is not altered. In Table S3 single point energy results at the M06/6-311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER and M062X/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER levels are reported. Even though the absolute potential energy barriers and reaction energies with the M06 and M06X functionals are somewhat different than the B3LYP values, the qualitative trends are maintained. In Table S4 the distances corresponding to the breaking and forming bonds of the HzproR transfers are given for the reactant, transition state structure and product for the H13proR abstraction for snapshots I, II, and IV, and for the reactant and product structures for snapshots II and V of the H16proR abstraction. As previously observed for the COX-2/AA28 and the COX-1/AA systems,42 the dispersion in the potential energy barriers for the H13proR transfer is mainly due to the dispersion in the geometries of the reactant structures because the transition state geometries are more similar. Those H-transfers are somewhat exothermic, although slightly endothermic for the H16proR abstraction in snapshot II, showing those reaction energies some correlation with the H-abstraction barriers. The optimized products are EPA radicals delocalized over the C11−C15 planar pentadienyl system for the H13proR transfer, and delocalized over the C14-C18 planar pentadienyl system for the H16proR transfer. In Table S5, the values for the two 17 dihedrals of each pentadienyl group are given. Also, we have included the C11-C12 and C14-C15 distances for the pentadienyl radical centered at C13, and the C14-C15 and C17-C18 distances for the pentadienyl radical centered at C16. The similarity of those two distances for each pentadienyl group demonstrates that both are delocalized radicals in all the snapshots, even though the corresponding dihedrals deviate a little from planarity. 3.3 Molecular dynamics simulation of the aspirin-acetylated COX-2/EPA complex For the aspirin-acetylated COX-2/EPA complex, we have only performed a MD simulation of 100 ns in the presence of the Tyr385-O· radical. The system reaches equilibrium during the trajectory according to the stable behaviour of the protein RMSD depicted in Figure S5. Again, the N-terminal domain is the most flexible part of the biomolecular complex, as can be seen in Figure S6. As indicated in the Introduction section, the main products of EPA biosynthesis by aspirin-acetylated COX-2 are 18R-HEPE and 18S-HEPE, which come from the H16-abstraction from EPA.16 However, we have also analyzed here the H13-abstraction process trying to reveal the molecular insights that could explain the regioselectivity of EPA catalysis by aspirin-acetylated COX-2. The evolution of the H13proX- OTyr385 and H16proX-OTyr385 distances has been recorded and plotted in Figure 5. The plot of the H13proR-OTyr385 distance shows an oscillatory behaviour taking an average value of 3.78 Å in the last 20 ns. The H13proS-OTyr385 distance also presents some oscillations at the beginning of the simulation, but then it becomes more stable around an average value of 3.40 Å. The H16proX hydrogens move further away after 40 ns due to a substrate rearrangement. The average H16proX-OTyr385 distances are larger than the average H13proX-OTyr385 distances most of the time. This result seems not to agree with the regioselectivity of the aspirin-acetylated COX-2 enzyme in EPA catalysis.16,17 In any case, further analysis is needed. All the calculated average distances for the MD simulation are collected in Table S6. Following, we filtered 10,000 frames of the 100 ns trajectory by taking one snapshot each 10 ps. According to the geometric conditions detailed in Section 3.1 to define the pre-catalytic structures for 18 the wild type COX-2/EPA complex, we then selected from those 10,000 frames a group of complexes ready to initiate the H13proX and/or the H16proX abstraction processes. The percentage of pre- catalytic structures found confirms that there are more structures suited for H13proX abstraction (35.0 % and 38.3 % for H13proR and H13proS transfers, respectively) than for the H16proX transfer (15.5 % for H16proR and 19.5 % for H16proS). Contrarily, the difference between the pre-catalytic structures for HZproR and HZproS transfers is not significant in the aspirin-acetylated COX-2/EPA complex. However, we decided to study the abstraction of the HZproS hydrogens to better compare with our previous study of the all-radical mechanism in the aspirin-acetylated COX-2/AA system.26 Those results demonstrated that the most favourable pathway for the formation of the main product (15R-HPETE) was initiated by H13proS abstraction. Figure 5. Evolution of the H13proX-OTyr385 and H16proX-OTyr385 distances along the 100ns MD simulation of the aspirin-acetylated COX-2/EPA complex in the presence of the Tyr385-O· radical. On the right side, the distances H16proX-OTyr385 are represented. On the left side, the H13proX-OTyr385 distances are shown. Next, we performed a clustering analysis based on the heavy-atom RMSD of the acetyl-Ser530 group, and five clusters were obtained. Then, an intersection was done between the pre-catalytic structures selected previously and those Ser530-based clusters. We will call snaphot I and snapshot II to the centroids of the intersection between the first and second cluster, respectively, and the 3830 pre- 19 catalytic structures chosen to initiate the H13proS transfer. The centroid corresponds to that structure with the smallest EPA heavy atoms RMSD with respect to the average position of those heavy atoms over all the structures of the intersection set. Likewise, snapshots V and VI will correspond to the centroids of the intersection between the first and second clusters, respectively, and the 1954 pre- catalytic structures corresponding to the abstraction of H16proS. Furthermore, two snapshots (that will be denominated snapshot III and IV) with H13proS near the oxygen acceptor of the Tyr385-O· radical have been selected randomly from the intersection of the pre-catalytic structures of H13proS with the first cluster of the aspirin-acetylated Ser530. Finally, two more snapshots (called snapshot VII and VIII) with H16proS near the Tyr385-O· radical have been selected randomly from the intersection of the pre- catalytic structures of H16proS with the first cluster of the aspirin-acetylated Ser530. The initial distances corresponding to the three atoms directly involved in the bond breaking and bond forming of the hydrogen abstractions (first step) of the all-radical mechanism for the eight selected snapshots are given in Table 4. It must be pointed out that in snapshots II and V, H16proS and H13proS, respectively, are also quite close to the tyrosyl radical. For that reason, both structures have also been clasified as suitable for the study of the H16proS and H13proS abstraction (see Table 4). Table 4. Initial distances (Å) corresponding to the three atoms directly involved in the hydrogen abstraction (first step) of the all-radical mechanism for the eight selected snapshots in the case of the aspirin-acetylated COX-2/EPA complex in the presence of the Tyr385-O· radical. H13proS H16proS Snapshots d(H-O) d(C-O) Snapshots d(H-O) d(C-O) I 2.82 3.81 V 2.86 3.49 II 2.92 3.85 II 3.31 2.86 III 2.73 3.33 VI 2.79 3.79 IV 2.60 3.25 VII 2.49 3.33 V 2.65 3.73 VIII 2.25 2.87 20 The orientation of the acetyl-Ser530 group in the first and second clusters of the aspirin-acetylated COX-2 system corresponds to the conformations of clusters II and IV of our previous study of the aspirin-acetylated COX-2/AA complex.26 In those clusters, the acetyl group of Ser530 adopts a parallel orientation with respect to the hydrophobic groove, and EPA is bound with an “L-shaped” conformation as in the wild type complex and also similarly to AA. In Figure 6, we have depicted structures I, IV, V, and VII. In the four snapshots, the carboxylate group of the EPA substrate interacts with Tyr355 and Arg120 (see Figure S7 and Table S2). Also, the structures ready for H13proS abstraction (snapshots I and IV) present a more extended binding mode than those corresponding to H16proS abstraction (snapshots V and VII). In Figure S8, an overlay between snapshots IV and VII shows the difference between those two binding modes. In comparison with the COX-2/EPA system, we observe a more “L-shaped” binding mode of EPA in the presence of the aspirin-acetylated Ser530. That is evident for the reactive structures of both H13proS and H16proS transfers, as shown in Figures S9 A and S9 B, respectively. The acetyl group of aspirin-acetylated Ser530 would block the disposition of the EPA ligand if it were bound as in the COX-2 active site. 21 Figure 6. Binding modes of the EPA substrate inside the hydrophobic groove of aspirin-acetylated COX-2 for snapshots I, IV, V and VII in the presence of the Tyr385-O· radical. 3.4 QM/MM calculations of the H13proS/H16proS abstraction processes in the aspirin-acetylated COX- 2/EPA complex According to the all-radical mechanism in the aspirin-acetylated COX-2/EPA system, we first studied the H13proS and H16proS abstraction transfers (first step of the all-radical mechanism) from EPA to the Tyr385-O· radical using QM/MM techniques. In the next section we will focus on the molecular oxygen addition to the formed EPA pentadienyl radical. For the sake of clarity, the reactions studied for the aspirin-acetylated COX-2/EPA system are schematized in Figure 7. Figure 7. Scheme of all the studied reaction pathways leading to 18R-HEPE and 18S-HEPE from EPA within aspirin-acetylated COX-2. Step 1 corresponds to the initial H-abstraction. Steps 11 and 12 correspond to the formation of 14-peroxy-eicosapentaenoic radical and 18-peroxy-eicosapentaenoic radical, respectively. The paths associated to high potential energy barriers are marked by a red cross. 22 For the hydrogen abstractions, first, we optimized the structures of the reactant complexes corresponding to the different snapshots. Next, we calculated the potential energy profiles of those H-abstraction reactions using a reaction coordinate defined as the distance corresponding to the forming bond (HzproS − OTyr385) for both the H13proS and H16proS atoms. The transition state and the product structures were taken from the scan and were also fully optimized. The potential energy barriers for the H13proS and H16proS abstraction reactions are given in Table 5 along with the reaction energies. Selected geometrical parameters of reactant, transition state structure and product geometries are given in Table S7. The dispersion on the potential energy barriers is correlated with the dispersion of the reactant geometries because the H-abstraction transition state structures present very similar structures for all the selected snapshots in accordance with previous studies on the same type of reaction.45,46 It is interesting to highlight that in the aspirin-acetylated COX-2/EPA complex, the H16proS abstraction energy barriers tend to be somewhat lower than the H13proS energy barriers. Also, the H16proS barriers in the aspirin-acetylated COX-2/EPA complex are clearly smaller than for the H16proR transfers in the wild type COX-2/EPA system. These results are in accordance with the formation of 18R-HEPE and 18S-HEPE as the main products from EPA in aspirin-acetylated COX-2.16,18 The reaction energies are slightly negative or positive except the one for snapshot II (H13proS) that is more endoergic, and that for snapshot VII that is the most exoergic. As seen in Table 5 the QM/MM single point energy results at the B3LYP/6-311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER level qualitatively agree with the fully optimized B3LYP/6-31G(d,p):AMBER ones. Also, the single point energy calculations reported in the last column of Table 5 show that the inclusion of dispersion corrections lowers the potential energy barriers43,44 even though the qualitative trends of those energy barriers are maintained. The reactant and transition state structures of the H16proS abstraction pathway for snapshot VII were relocalized at the B3LYP/6-31G(d,p)-GD3BJ:AMBER level to check the effect of including dispersion corrections in the optimization procedure. The recalculated potential energy barrier has a value of 22.3 kcal/mol, which is only 1 kcal/mol lower than the B3LYP/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER energy barrier including a posteriori GD3BJ 23 dispersion corrections. Therefore, the computational strategy of including dispersion corrections via single point energy calculations with a large basis set gives similar results for the H-abstraction barrier than optimizing with dispersion corrections but with a smaller basis set. In Table S8 single point energy results at the M06/6-311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER and M062X/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER levels are included. Even though the absolute potential energy barriers and reaction energies with the M06 and M06X functionals are somewhat different than the B3LYP values, the qualitative trends are not altered. In Table S9, the values for the two dihedrals associated with the radicals centered at C13 or C16 are given for the corresponding products of the abstraction reaction. Those radicals are nearly planar in all the cases except for snaphots II and V corresponding to the H13proS abstraction. This means that the products of the reaction profiles of snapshots I, III and IV for H13proS, and all those corresponding to H16proS are structures stabilized by conjugation in a pi system of five electrons delocalized over the corresponding pentadienyl groups. Contrarily, the products of snapshots II and V (H13proS) should be better described as allyl radicals only delocalized over C13-C15. The pentadienyl radicals for H13proS present very similar C11-C12/C14-C15 bond distances and also the C14-C15/C17-C18 bond distances are very close for H16proS product structures. The allyl radicals for snapshots II and V (H13proS) have shorter C11- C12 bond distances. Note that the most endoergic H-transfer corresponds to the formation of the allyl radical of snapshot II (H13proS) whereas the most exoergic corresponds to the most planar pentadienyl radical of snapshot VII (H16proS). It is worth mentioning that the disposition of the pentadienyl group of the EPA radicals in the aspirin-acetylated COX-2 active site (See Figure 8 for snapshots I and VII), after either H13proS or H16proS abstractions, is very similar to that of the corresponding radicals for AA. In our previous paper on the aspirin-acetylated COX-2/AA system,26 those pentadienyl radicals successfully led to the formation of 15R-HETE after oxygen addition on C13 in agreement with experiments.47,48 24 Table 5. Potential energy barriers and reaction energies (in kcal/mol) at the B3LYP/6-31G(d,p):AMBER level corresponding to the H13proS and H16proS hydrogen abstractions for those snapshots that have been able to reach the products in the aspirin-acetylated COX-2/EPA complex, as well as the corresponding values obtained by means of QM/MM single point energy calculations at the B3LYP/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER level, and including a posteriori GD3BJ dispersion corrections. B3LYP/6-311+G B3LYP/6-311+G Snapshots B3LYP/6-31G(d,p) (2d,2p) (2d,2p) - GD3BJ H13proS ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E I 27.7 -0.9 29.1 -1.8 26.6 -2.6 II 27.5 9.7 28.2 8.8 27.3 6.4 III 27.4 1.6 28.6 1.5 27.3 -2.4 IV 27.0 2.5 27.9 2.1 23.4 -1.6 V 24.3 5.2 25.0 4.4 25.4 4.9 H16proS ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E V 20.4 -3.8 22.7 -3.9 16.6 -7.4 II 26.6 1.8 27.5 1.7 26.0 -0.9 VI 26.9 2.9 29.2 2.3 24.5 -1.7 VII 24.6 -12.0 26.1 -12.3 23.3 -13.8 VIII 25.7 2.8 26.7 2.1 22.2 -2.4 Figure 8. Overlay between the EPA radicals formed after H13proS and H16proS abstractions in snapshots I (beige) and VII (gray), respectively, in the aspirin-acetylated COX-2 active site. Only hydrogens of EPA at C11-C15 and C14-C18 are depicted. 25 After the corresponding H-abstraction, the formation of the delocalized C11−C15 and C14-C18 planar pentadienyl radicals would allow the O2 attack on C11 and C15, or on C14 and C18. As our main objective was the study of the viability of 18R/S-HEPE formation by aspirin-acetylated COX-2, we decided to investigate the second step of the all-radical mechanism once H16proS is abstracted, and oxygen molecules enter the active site. Thus, snapshots VII and VIII were selected for the calculation of O2 addition energy profiles on C14 or C18. 3.5 QM/MM calculations of the oxygen addition reaction in the aspirin-acetylated COX-2/EPA complex The second step of the all-radical mechanism in the COX-2/EPA complex consists of the addition of one oxygen molecule to the C14-C18 EPA pentadienyl radical formed once H16proS is abstracted. As mentioned above, the EPA radicals for all the H16proS snapshots present a quite planar pentadienyl fragment so corresponding to delocalized radicals over the five carbon atoms from C14 to C18. Thus, the oxygen molecule could either attack to C14 or C18 because the terminal carbon atoms present similar spin densities that are both bigger than for the central carbon. We have explored different initial locations for the O2 molecule around C14 or C18 at the conformation of the EPA pentadienyl radicals for snapshot VII and VIII, respectively (see Figure 9). The O2 molecules around C14 or C18, taken as origin of coordinates, have been initially placed along the x, y and z Cartesian axes and along the bisector axes contained in the xy, xz and yz planes. 18 O2 molecules have been set to distances of 3.0 and 3.5 Å. 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 or C18 (for the sake of example, the 8 starting points for the case of the O2 addition to C18 in snapshot VII are depicted in 26 Figure S10). The reaction coordinate was taken as the distance from the attacking oxygen of the O2 molecule to C14 or C18 of the EPA pentadienyl radical. The potential energy barriers for the O2 addition reaction coordinates that were succesful leading to the peroxyl radical products are collected in Tables 6 and 7 for the addition to C14 or C18, respectively. Figure 9. Conformations of the delocalized C14–C18 pentadienyl radicals formed inside the hydrophobic groove of aspirin-acetylated COX-2 after hydrogen abstraction of H16proS. Snapshot VII (up) is represented in blue, and snapshot VIII (bottom) in green. Table 6. Potential energy barriers and reaction energies (in kcal/mol) corresponding to the O2 addition to C14 at the B3LYP/6-31G(d,p):AMBER level for all the pathways that have been able to reach the products for snapshot VII in the aspirin-acetylated COX-2/EPA complex, as well as the corresponding values obtained by means of QM/MM single point energy calculations at the B3LYP/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER level, including GD3BJ dispersion corrections added 27 afterwards. Antara and supra are the shortening word for antarafacial and suprafacial attack, respectively. B3LYP/6-311+G B3LYP/6-311+G B3LYP/6-31G(d,p) (2d,2p) (2d,2p) - GD3BJ Stereochemistry O2 attack ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E C14 – Snap. VII Antara R 16.4 -10.9 15.7 -9.2 18.8 -11.5 Antara R 11.4 -2.0 12.6 -0.2 12.8 1.9 Antara R 11.3 -2.4 12.4 -0.5 12.1 1.5 Antara R 22.7 -5.0 27.4 5.0 25.8 4.6 Antara R 11.3 -2.2 12.5 -0.4 12.6 1.5 Supra S 15.4 3.3 15.8 4.3 9.1 -0.6 First, it is interesting to note that the O2 attacks can take place following a suprafacial or an antarafacial pathway. In a suprafacial attack, the O2 molecule approaches the pentadienyl radical by the same side of the Tyr385-O· radical, whereas in an antarafacial approach the O2 molecule approaches the pentadienyl radical by the opposite side of the Tyr385-O· radical. For the addition to C14, five of the six addition pathways are antarafacial and only one suprafacial. The stereochemistry of C14 at the final peroxyl radicals formed is correlated with the antarafacial or suprafacial character of the O2 addition. So, when the addition is antarafacial, C14 presents an R stereochemistry, but if the O2 attack is suprafacial C14 has an S stereochemistry. As for the potential energy barriers for O2 addition at C14, all of them are too high in comparison with the energy barriers calculated for the same reaction step in COX-2/AA and aspirin-acetylated COX-2/AA systems that we studied previously.26,28 In Table S10 the energy barriers obtained by means of QM/MM single point energy calculations at the M06/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER and M062X/6-311+G(2d,2p):AMBER//B3LYP/6- 31G(d,p):AMBER levels are reported. Those results confirm that the oxygen addition reactions at C14 are too high for being competitive. In Table 7, the results for the O2 additions to C18 are collected for snapshots VII and VIII. In this case, when the addition is antarafacial C18 presents an S stereochemistry, but if the O2 attack is suprafacial 28 C18 has an R stereochemistry. For the non suprafacial/antarafacial approaches (see Figure 10 as example), the final peroxyl radical at C18 is R. The potential energy barriers for both the antarafacial and the suprafacial attacks are quite lower than in the case of the O2 addition to C14. Therefore, formation of the peroxyl radical intermediates leading to both 18S-HEPE and 18R-HEPE via this mechanism in aspirin-acetylated COX-2 appears to be possible. In any case, the suprafacial attacks (giving a final R stereochemistry at the peroxyl product) present in general lower barriers. As for the hydrogen abstractions, in the case of the oxygen additions the QM/MM single point energy calculations lead qualitatively to the same conclusions as with the fully optimized B3LYP/6- 31G(d,p):AMBER values. In Table S11 the energy barriers obtained by means of QM/MM single point energy calculations at the M06/6-311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER and M062X/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER levels are reported. Those results confirm that the oxygen addition reactions at C18 are viable, showing lower energy barriers than additions at C14. Note that the stereochemistry of the final peroxyl radical can depend on each snapshot. This is because the delocalized C14−C18 pentadienyl radical adopts, after the H-abstraction, a different conformation in each snapshot. If the pentadienyl group presents its S face to the antarafacial side, an O2 antarafacial attack will lead to an S stereochemistry at the peroxyl radical. The R stereochemistry will result from the suprafacial O2 attack to the R face of the pentadienyl group. It is worth mentioning that when the oxygen attack is suprafacial giving an R stereochemistry at C18, the peroxo group in the radical intermediates faces the OH group of Tyr385, even forming in some cases a hydrogen bond between an oxygen atom of the peroxo group at C18 and the OH group of Tyr385. Conversely, when the oxygen attack is antarafacial giving an S stereochemistry at C18, the peroxo group in the radical intermediate points to the opposite side of the OH group of Tyr385 (see Figure 11). Therefore, we can see that the 18R-peroxyl radicals coming from the suprafacial oxygen addition are already well prepared for the back-hydrogen transfer from Tyr385 to the terminal oxygen of the peroxo group, to form the 18R-hidroperoxy-eicosapentaenoic (18R-HpEPE). On the contrary, 29 the 18S-peroxyl radicals coming from the antarafacial oxygen addition would still require a previous rotation of the peroxo group that would involve a significant potential energy barrier to form 18S- HpEPE. All that explains why, after reduction of the hydroperoxides, the predominant 18-HEPE isomer is 18R-HEPE. Table 7. Potential energy barriers and reaction energies (in kcal/mol) corresponding to the O2 addition to C18 at the B3LYP/6-31G(d,p):AMBER level for all the pathways that have been able to reach the products for snapshots VII and VIII in the aspirin-acetylated COX-2/EPA complex, as well as the corresponding values obtained by means of QM/MM single point energy calculations at the B3LYP/6- 311+G(2d,2p):AMBER//B3LYP/6-31G(d,p):AMBER level, including GD3BJ dispersion corrections added afterwards. Antara and supra are the shortening word for antarafacial and suprafacial attack, respectively. B3LYP/6-311+G B3LYP/6-311+G B3LYP/6-31G(d,p) (2d,2p) (2d,2p) - GD3BJ Stereochemistry O2 attack ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E C18 –Snap. VII Antara S 6.7 -18.7 6.5 -18.8 3.8 -15.5 Non R 4.1 -11.5 3.8 -11.4 5.3 -9.0 Supra/Antara Non R 4.2 -11.6 3.8 -11.5 5.2 -9.1 Supra/Antara Non R 4.0 -11.7 3.7 -11.5 5.2 -9.1 Supra/Antara Supra R 1.0 -11.0 1.2 -10.4 0.2 -9.4 Stereochemistry O2 attack ∆E‡ ∆E ∆E‡ ∆E ∆E‡ ∆E C18 –Snap. VIII Antara S 6.9 -5.7 5.9 -2.9 2.5 -6.3 Antara S 0.7 -14.2 0.7 -12.7 0.1 -11.6 Supra R 2.7 0.0 2.5 -7.3 2.4 -6.9 Supra R 2.1 -14.6 7.2 -14.2 8.1 -12.7 30 Figure 10. Example of Non Suprafacial/Antarafacial O2 addition to C18 for snapshot VII in the aspirin- acetylated COX-2/EPA complex. 31 Figure 11. Formation of the peroxyl radical intermediates then leading to 18R-HEPE and 18S-HEPE in the aspirin-acetylated COX-2/EPA complex. Up: 18R-peroxyl radicals from the suprafacial oxygen addition to snapshot VII (blue) or snapshot VIII (green). Bottom: 18S-peroxyl radicals from the antarafacial oxygen addition to snapshot VIII (green). 4. Conclusions Acute inflammation is a very important defense mechanism of host against a variety of harmful stimuli, that consists of an onset phase, followed by a resolution phase. This resolution phase is regulated by a number of specialized pro-resolving lipid mediators (SPMs) produced by human cells called macrophages. If resolution does not work well, the acute inflammation becomes a chronic inflammation, that is linked to a very wide variety of human diseases, including Coronavirus Disease 2019 (COVID-19). A number of nonsteroidal anti-inflammatory drugs (NSAIDs), like aspirin, ibuprofen or diclofenac, have been found to be therapeutic drugs for the onset phase of the inflammation. However, most of these traditional drugs, although diminishing the onset phase by inhibiting COX-2, block the biosynthesis of SPMs as well, hence delaying resolution. In this sense, aspirin is a very important exception. Thus, aspirin-acetylated COX-2 triggers the biosynthesis of several SPMs from some polyunsaturated fatty acids. In particular, aspirin-acetylated COX-2 rather than COX-2, is able to initiate the biosynthesis of E-series resolvins from EPA, by means of the initial production of 18-hydroperoxy-eicosapentaenoic acid (18-HpEPE), that is subsequently reduced by a peroxidase to 18-HEPE. In this paper we have combined molecular dynamics simulations and QM/MM calculations to understand the molecular details of this very intriguing behaviour. We have shown that in the COX-2/EPA complex EPA presents a “bulged L-shaped” binding mode. The potential energy barriers for the H13 abstraction are clearly lower than for the H16 abstraction. It is worth pointing out that the H13-abstraction energy barriers for EPA are, on average, higher than those 32 for AA, in agreement with the experimental results that have proved that EPA is a poorer substrate of COX-2 than AA. Since the H16 abstraction is the one required for the formation of 18-HpEPE, our results indicate that this product should be clearly marginal because of the action of COX-2 on EPA. In comparison with the COX-2/EPA system, EPA adopts a more “L-shaped” binding mode in the presence of the aspirin-acetylated Ser530. The acetyl group of aspirin-acetylated Ser530 would block the disposition of the EPA ligand if it were bound as in the COX-2 active site. In the aspirin-acetylated COX-2/EPA complex, the H16proS abstraction energy barriers tend to be somewhat lower than the H13proS energy barriers. Also, the H16proS barriers in the aspirin-acetylated COX-2/EPA complex are clearly smaller than for the H16 transfers in the wild type COX-2/EPA system. In the following step, the O2 addition to C18 turns out to be very favoured versus the corresponding addition to C14. These results are in accordance with the formation of 18R-HEPE and 18S-HEPE as the main products from EPA in aspirin-acetylated COX-2. When the oxygen attack is suprafacial giving an R stereochemistry at C18 with somewhat lower energy barriers, the peroxo group in the radical intermediates faces the OH group of Tyr385. Conversely, when the oxygen attack is antarafacial giving an S stereochemistry at C18, the peroxo group in the radical intermediate points to the opposite side of the OH group of Tyr385. Therefore, the 18R-peroxyl radicals coming from the suprafacial oxygen addition are already well prepared for the back-hydrogen transfer from Tyr385 to the terminal oxygen of the peroxo group, to form 18R-HpEPE. All that explains why, after reduction of the hydroperoxides, the predominant 18- HEPE isomer is 18R-HEPE, although 18S-HEPE is also produced. These two isomers are the precursors of the E-series resolvins. Electronic Supplementary Information RMSDs corresponding to the MD simulations of the COX-2/EPA Michaelis complex; Binding modes of the EPA substrate inside the hydrophobic groove of COX-2; Potential energy profiles for H abstractions in the COX-2/EPA complex; RMSDs corresponding to the MD simulation of the aspirin-acetylated COX- 2/EPA Michaelis complex; Binding modes of the EPA substrate inside the hydrophobic groove of 33 aspirin-acetylated COX-2; Overlays between snapshots of COX-2/EPA and aspirin-acetylated COX- 2/EPA complexes; Starting points for the O2 addition to C18 in snapshot VII in the aspirin-acetylated COX-2/EPA complex; Average distances along the MD simulations of the COX-2/EPA Michaelis complex; Main hydrogen bonds at the entrance of the hydrophobic groove for representative pre- catalytic snapshots of COX-2/EPA and aspirin-acetylated COX-2/EPA complexes; Interatomic distances directly involved in the H13proR and H16proR abstractions and potential energy barriers at the B3LYP/6- 31G(d,p):AMBER optimized reactant, transition state structure and product corresponding to the selected snapshots of the COX-2/EPA Michaelis complex; Distances and dihedral angles of the optimized radicals produced as products of the first H-abstraction step for snapshots I, II, IV, and V of the COX-2/EPA system; Average distances along the MD simulation of the aspirin-acetylated COX- 2/EPA Michaelis complex; Interatomic distances directly involved in the H13proS and H16proS abstractions and potential energy barriers at the optimized reactant, transition state structure and product corresponding to the selected snapshots of the aspirin-acetylated COX-2/EPA Michaelis complex; Distances and dihedral angles of the optimized radicals produced as products of the first H-abstraction step for the selected snapshots of the aspirin-acetylated COX-2/EPA system. Author contributions P. A. F., A.G-L and J.M.L. contributed to the conceptualization of the research goals and aims; designed the methodology; supervised the overall research activity and its execution and drafted the original manuscript as well as the review and editing of the revised version. A.C-P. and A.P. carried out the investigation and the formal analysis; A.C-P. prepared the visualization/data presentation. Acknowledgements The authors acknowledge the Spanish “Ministerio de Ciencia, Innovación y Universidades” for Grant CTQ2017-83745-P and the Fundação para a Ciência e a Tecnologia for financial support through the grant PTDC/QUI-OUT/1401/2020. 34 References 1 D. Panigrahy, M. M. Gilligan, S. Huang, A. Gartung, I. Cortés-Puch, P. J. Sime, R. P. Phipps, C. N. Serhan and B. D. Hammock, Cancer and Metastasis Reviews, 2020, 39, 337–340. 2 U. N. Das, Archives of Medical Research, 2020, 51, 282–286. 3 C. N. Serhan and N. A. Petasis, Chemical Reviews, 2011, 111, 5922–5943. 4 C. N. Serhan, N. Chiang and T. E. van Dyke, Nature Reviews Immunology, 2008, 8, 349–361. 5 C. N. Serhan, N. Chiang and J. Dalli, Molecular Aspects of Medicine, 2018, 64, 1–17. 6 C. N. Serhan and N. Chiang, Current Opinion in Pharmacology, 2013, 13, 632–640. 7 C. N. Serhan, Nature, 2014, 510, 92–101. 8 M. Arita, F. Bianchini, J. Aliberti, A. Sher, N. Chiang, S. Hong, R. Yang, N. A. Petasis and C. N. Serhan, Journal of Experimental Medicine, 2005, 201, 713–722. 9 L. G. Cleland, M. J. James and S. M. Proudman, Arthritis research & therapy, 2006, 8, 202. 10 L. G. Cleland, M. J. James and S. M. Proudman, Arthritis Research and Therapy, 2006, 8, 2006. 11 M. Wada, C. J. DeLong, Y. H. Hong, C. J. Rieke, I. Song, R. S. Sidhu, C. Yuan, M. Warnock, A. H. Schmaier, C. Yokoyama, E. M. Smyth, S. J. Wilson, G. A. FitzGerald, R. M. Garavito, X. S. De, J. W. Regan and W. L. Smith, Journal of Biological Chemistry, 2007, 282, 22254–22266. 12 N. P. Sharma, L. Dong, C. Yuan, K. R. Noon and W. L. Smith, Molecular Pharmacology, 2010, 77, 979–986. 13 C. N. Serhan, C. B. Clish, J. Brannon, S. P. Colgan, N. Chiang and K. Gronert, Journal of Experimental Medicine, 2000, 192, 1197–1204. 14 E. Tjonahen, S. F. Oh, J. Siegelman, S. Elangovan, K. B. Percarpio, S. Hong, M. Arita and C. N. Serhan, Chemistry and Biology, 2006, 13, 1193–1202. 15 C. N. Serhan, S. Hong, K. Gronert, S. P. Colgan, P. R. Devchand, G. Mirick and R. L. Moussignac, Journal of Experimental Medicine, 2002, 196, 1025–1037. 16 S. F. Oh, T. W. Vickery and C. N. Serhan, Biochimica et Biophysica Acta - Molecular and Cell Biology of Lipids, 2011, 1811, 737–747. 17 S. F. Oh, P. S. Pillai, A. Recchiuti, R. Yang and C. N. Serhan, Journal of Clinical Investigation, 2011, 121, 569–581. 18 C. N. Serhan, FASEB Journal, 2017, 31, 1273–1288. 19 A. J. Vecchio, D. M. Simmons and M. G. Malkowski, Journal of Biological Chemistry, 2010, 285, 22152–22163. 20 J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser and C. Simmerling, Journal of Chemical Theory and Computation, 2015, 11, 3696–3713. 21 D. A. Case, R. M. Betz, D. S. Cerutti, T. E. Cheatham III, T. A. Darden, R. E. Duke, T. J. Giese, H. Gohlke, A. W. Goetz, N. Homeyer, S. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T. S. Lee, S. 35 LeGrand, P. Li, C.Lin, T. Luchko, R. Luo, B. Madej, D. Mermelstein, K. M. Merz, G. Monard, H. Nguyen, H. T. Nguyen, I. Omelyan, A. Onufriev, D. R. Roe, A. Roitberg, C. Sagui, C. L. Simmerling, W. M. Botello-Smith, J. Swails, R. C. Walker, J. Wang, R. M. Wolf, X. Wu, L. Xiao and P. A. Kollman, (2016) AMBER 2016 University of California, San Francisco . 22 K. Shahrokh, A. Orendt, G. S. Yost and T. E. Cheatham, Journal of Computational Chemistry, 2012, 33, 119–133. 23 C. I. Bayly, P. Cieplak, W. D. Cornell and P. A. Kollman, Journal of Physical Chemistry, 1993, 97, 10269–10280. 24 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, Journal of Computational Chemistry, 2004, 25, 1157–1174. 25 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, The Journal of Chemical Physics, 1983, 79, 926–935. 26 A. Cebrián-Prats, À. González-Lafont and J. M. Lluch, ACS Catalysis, 2020, 10, 138–153. 27 M. J. Lucido, B. J. Orlando, A. J. Vecchio and M. G. Malkowski, Biochemistry, 2016, 55, 1226– 1238. 28 A. Cebrián-Prats, À. González-Lafont and J. M. Lluch, ACS Omega, 2019, 4, 2063–2074. 29 M. Deserno and C. Holm, Journal of Chemical Physics, 1998, 109, 7678–7693. 30 H. Lei, N. A. Baker and X. Li, Proceedings of the National Academy of Sciences of the United States of America, 2016, 113, 14183–14188. 31 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola and J. R. Haak, The Journal of Chemical Physics, 1984, 81, 3684–3690. 32 S. le Grand, A. W. Götz and R. C. Walker, Computer Physics Communications, 2013, 184, 374– 380. 33 R. Salomon-Ferrer, A. W. Götz, D. Poole, S. le Grand and R. C. Walker, Journal of Chemical Theory and Computation, 2013, 9, 3878–3888. 34 F. Maseras and K. Morokuma, Journal of Computational Chemistry, 1995, 16, 1170–1179. 35 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. a. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. a. Petersson, H. Nakatsuji, X. Li, M. Caricato, a. v. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. v. Ortiz, a. F. Izmaylov, J. L. Sonnenberg, Williams, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. a. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. a. Keith, R. Kobayashi, J. Normand, K. Raghavachari, a. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian, Inc., Wallingford CT, 2009. 36 A. D. Becke, The Journal of Chemical Physics, 1993, 98, 5648–5652. 37 C. Lee, W. Yang and R. G. Parr, Physical Review B, 1988, 37, 785–789. 36 38 S. H. Vosko, L. Wilk and M. Nusair, Canadian Journal of Physics, 1980, 58, 1200–1211. 39 W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, Journal of the American Chemical Society, 1996, 118, 2309–2309. 40 A. R. Calixto, M. J. Ramos and P. A. Fernandes, Journal of Chemical Theory and Computation, 2017, 13, 5486–5495. 41 S. Grimme, S. Ehrlich and L. Goerigk, Journal of Computational Chemistry, 2011, 32, 1456– 1465. 42 C. Z. Christov, A. Lodola, T. G. Karabencheva-Christova, S. Wan, P. v. Coveney and A. J. Mulholland, Biophysical Journal, 2013, 104, L5–L7. 43 R. Lonsdale, J. N. Harvey and A. J. Mulholland, Journal of Chemical Theory and Computation, 2012, 8, 4637–4645. 44 R. Lonsdale, J. N. Harvey and A. J. Mulholland, The Journal of Physical Chemistry Letters, 2010, 1, 3232–3237. 45 P. Saura, R. Suardíaz, L. Masgrau, J. M. Lluch and À. González-Lafont, ACS Catalysis, 2014, 4, 4351–4363. 46 P. Saura, I. Kaganer, D. Heydeck, J. M. Lluch, H. Kühn and À. González-Lafont, Chemistry - A European Journal, 2018, 24, 962–973. 47 G. Xiao, A. L. Tsai, G. Palmer, W. C. Boyar, P. J. Marshall and R. J. Kulmacz, Biochemistry, 1997, 36, 1836–1845. 48 J. A. Gimenez-Bastida, W. E. Boeglin, O. Boutaud, M. G. Malkowski and C. Schneider, FASEB Journal, 2019, 33, 1033–1041. 37