RSC Advances This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online PAPER View Journal | View Issue Unraveling how the Gly526Ser mutation arrests prostaglandin formation from arachidonic acid Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. Cite this: RSC Adv., 2020, 10, 986 catalyzed by cyclooxygenase-2: a combined molecular dynamics and QM/MM study† Adrián Suñer-Rubio, a Anna Cebrián-Prats, a Àngels González-Lafont ab and José M. Lluch *ab Cyclooxygenases (COXs) are the enzymes responsible for the biosynthesis of prostaglandins, eicosanoids that play a major role in many physiological processes. Particularly, prostaglandins are known to trigger inflammation, and COX-2, the enzyme isoform associated with this inflammatory response, catalyzes the cyclooxidation of arachidonic acid, leading to prostaglandin G2. For this reason, COX-2 has been a very important pharmacological target for several decades now. The catalytic mechanism of COX-2, a so- called all-radical mechanism, consists of six chemical steps. One of the most intriguing aspects of this mechanism is how COX-2 manages to control the regio- and stereospecificity of the products formed at each step. Mutagenesis experiments have previously been performed in an attempt to find those hot- spot residues that make such control possible. In this context, it is worth mentioning that in experiments Received 28th October 2019 Accepted 12th December 2019 with the Gly526Ser COX-2 mutant, prostaglandins were not detected. In this paper, we have combined molecular dynamics simulations and quantum mechanics/molecular mechanics calculations to analyze DOI: 10.1039/c9ra08860a how the COX-2 catalytic mechanism is modified in the Gly526Ser mutant. Therefore, this study provides rsc.li/rsc-advances new insights into the COX-2 catalytic function. (COX-2), which leads to the formation of the eicosanoid pros- 1 Introduction taglandins. Prostaglandins control many pathological events, Lipid mediators1 are bioactive lipids involved in many physio- such as inammation, fever, pain, cardiovascular diseases, and logical processes, and are generated through specic biosyn- cancer.1 thetic pathways. They bind to their corresponding G protein- AA is released from membrane glycerophospholipids. COX-2 coupled receptors and transmit signals to target cells. The is a membrane-associated homodimeric bifunctional hemo- synthesis of many lipid mediators begins with the oxygenation protein.3 Only one monomer acts as a catalyst at a given time.4 of polyunsaturated fatty acids by some enzymes.2 One of the The determination of the corresponding crystal structure3 has most well known and important cases is the dioxygenation of shown that in the murine cyclooxygenase active site of the arachidonic acid (AA, 20:4 n-6) by the enzyme cyclooxygenase-2 catalytic monomer the AA carboxylate lies near the side chains a Departament de Quı́mica, Universitat Autònoma de Barcelona, 08193 Bellaterra, populated clusters; evolution of the C13 pro-S hydrogen–O(Tyr385) and Barcelona, Spain. E-mail: JoseMaria.Lluch@uab.cat C13–O(Tyr385) distances during 250 ns of the last MD simulation of the b Institut de Biotecnologia i de Biomedicina (IBB), Universitat Autònoma de Gly526Ser COX-2/AA Michaelis complex; evolution of the C13 pro-S Barcelona, 08193 Bellaterra, Barcelona, Spain hydrogen–O(Tyr385) and C13–O(Tyr385) distances during a new 100 ns MD † Electronic supplementary information (ESI) available: The QM region of each simulation of the Gly526Ser COX-2/AA Michaelis complex; distances between particular reaction step; evolution of the C13 pro-S hydrogen–O(Tyr385) distance side chain carbon atoms of selected residues of snapshot 1 extracted from the during 100 ns of the rst MD simulation of the Michaelis complex Gly526Ser MD simulation of the COX-2/AA Michaelis complex or the centroid snapshot COX-2/AA; evolution of the C13–O(Tyr385) distance during 100 ns of the rst of the most populated cluster of the last MD simulation of the Gly526Ser MD simulation of the Michaelis complex Gly526Ser COX-2/AA; distance COX-2/AA Michaelis complex, and arachidonic acid bound in the between the hydrogen of the Ser526 OH group and the oxygen atom of the corresponding active site active site; distances corresponding to the three backbone of Met522 during 10 ns (equilibration) + 20 ns (production) MD atoms that directly participate in the breaking/forming of bonds for the simulation of the Gly526Ser COX-2 mutated enzyme without AA; comparison reactant, transition state structure and product corresponding to the C13 between the structure corresponding to snapshot 1 extracted from the MD pro-S hydrogen abstraction by the tyrosyl radical; evolution of the key simulation of the COX-2/AA Michaelis complex and the centroid snapshot of distances during C13 pro-S hydrogen abstraction by the tyrosyl radical; the most populated cluster of the last MD simulation of the Gly526Ser structure obtained aer the 8,12-cyclization to give a bicyclo endoperoxide COX-2/AA Michaelis complex; comparison between the centroid of the most which occurs in wild-type COX-2. See DOI: 10.1039/c9ra08860a 986 | RSC Adv., 2020, 10, 986–997 This journal is © The Royal Society of Chemistry 2020 View Article Online Paper RSC Advances of Arg120 and Tyr355, whereas the u-end is positioned in the hydrophobic groove above Ser530, exhibiting an “L-shaped” binding conformation. A tyrosyl radical is generated as a consequence of one-electron transfer from Tyr385 to the oxy- This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. ferryl protoporphyrin IX radical cation formed in the peroxidase active site aer activation by an alkyl hydroperoxide of its heme group. The C13 carbon atom of AA is positioned 2.95 Å below the phenolic oxygen of Tyr385. The nascent tyrosyl radical starts the cyclooxygenase reaction by abstraction of the C13 pro-S hydrogen atom when AA occupies the cyclooxygenase active Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. site.5,6 Aerwards, addition of two oxygen molecules converts AA to prostaglandin G2 (PGG2). Later, PGG2 is released to the peroxidase active site, where its 15-hydroperoxyl group is reduced to give prostaglandin H2 (PGH2). The biosynthesis of PGG2 in the cyclooxygenase active site is achieved by a six-step all-radical mechanism2,7–13 (see Scheme 1) Fig. 1 Evolution of the C13 pro-S hydrogen–O(Tyr385) (red line) and derived from proposals by Hamberg and Samuelsson. It C13–O(Tyr385) (blue line) distances over 100 ns of the first MD simu- successively consists of: a C13 pro-S hydrogen abstraction from lation of the Gly526Ser COX-2/AA Michaelis complex. AA by the tyrosyl radical 385, forming a planar delocalized C11– C15 pentadienyl radical (step 1); an antarafacial molecular oxygen addition at C11 to give an 11R peroxyl radical (step 2); atoms). It is simultaneously a exible molecule (owing to its a 9,11-cyclization leading to a C8-radical cyclic endoperoxide multiple single bonds) and a rigid molecule (owing to its four (step 3); an 8,12-cyclization yielding a bicyclo endoperoxide and C]C double bonds). The role of many active site residues of COX- a C13–C15 allyl radical (step 4); a second antarafacial molecular 2 is to make this control possible. A number of mutagenesis oxygen addition at the 15S position (step 5); and a nal back experiments have been performed to better understand the way transfer of hydrogen from Tyr385 to the peroxyl radical at C15 they act.14–16 Particularly interesting is the case of the Gly526Ser (step 6), leading to the formation of PGG2. mutant human COX-2.16 This mutation does not appear to inter- One of the most amazing features of the six-step all-radical fere with the rst three reaction steps. However, no formation of mechanism is the extremely efficient way how COX-2 manages prostaglandins was detected, 11R-hydroperoxyeicosatetraenoic to control the stereochemistry and regiochemistry of the reactions acid (11R-HPETE) becoming the main product. Brash and co- it catalyzes. AA is a very long substrate (it contains 20 carbon workers16 proposed that the 8,12-cyclization in the Gly526Ser Scheme 1 The all-radical mechanism for the formation of PGG2 from arachidonic acid catalyzed by wild-type COX-2. This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 986–997 | 987 View Article Online RSC Advances Paper This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. Fig. 2 A typical structure corresponding to the second half of the first MD simulation of the Gly526Ser COX-2/AA Michaelis complex. mutant is sterically hindered by misalignment of the substrate the protein surface to neutralize the total charge of the system. carbon chains. So, instead of prostaglandins they found novel An almost cubic box (100.2 Å  113.4 Å  92.5 Å) of TIP3P water products, which were identied as 8,9-11,12-diepoxy-13R- (or 15R)- molecules20 was constructed, with a minimum cutoff distance hydroperoxy derivatives of arachidonic acid. of 15 Å between the atoms of the enzyme and the edge of the The aim of this paper is to provide a detailed molecular view of periodic box. All the water molecules closer than 2.2 Å to any the way how an apparently little signicant mutation is able to atom of the enzyme or AA were removed. In all, the system fully arrest the PGG2 biosynthesis, instead yielding diepoxy contained about 80 000 atoms. products. Thus, we have combined molecular dynamics (MD) simulations with quantum mechanics/molecular mechanics 2.2 Parameterization (QM/MM) calculations to analyze how the Gly526Ser mutant of First of all, the parameters of the heme group were taken from COX-2 alters the six-step all-radical mechanism shown in Scheme a recent study.21 However, the corresponding partial charges of 1 that holds for wild-type COX-2, so contributing to the under- the ferric penta-coordinated high spin group (Fe3+) were standing of the very important catalytic function of COX-2. calculated using the Restrained Electrostatic Potential (RESP)22 method at the HF/6-31G(d,p) level of theory. In addition, the 2 Computational methods 2.1 System setup In a previous paper we studied13 the mechanism that governs the biosynthesis of PGG2 catalyzed by wild-type COX-2. That system was built from the murine COX-2/AA complex (PDB code 3HS5)3, replacing the Co3+ complexed to the protoporphyrin IX group by Fe3+. The heme prosthetic group was tethered to the axial histidine (His388). Only the monomer B (the one where AA adopts an “L-shaped” productive binding conguration, with the C13 carbon atom of AA near Tyr385) was modeled in the cyclooxygenase active site. Aer a 100 ns MD simulation a pair of snapshots was selected, and we showed that both were able to follow the all-radical mechanism displayed in Scheme 1, even- tually yielding PGG2. In this paper, we have taken one of those snapshots (snapshot 1), and we have replaced the Gly526 residue by a serine using the CHIMERA17 program. All the amino acid residues of the protein but His388 and Tyr385 were described using the AMBER ff14SB force eld.18 All water molecules were deleted and the system was solvated again by Fig. 3 Evolution of the C13 pro-S hydrogen–O(Tyr385) (red line) and means of the AMBER18 19 tleap program. A pH of 7 was C13–O(Tyr385) (blue line) distances over 100 ns of the last MD simu- assumed for the titratable residues. Five Na+ ions were added at lation of the Gly526Ser COX-2/AA Michaelis complex. 988 | RSC Adv., 2020, 10, 986–997 This journal is © The Royal Society of Chemistry 2020 View Article Online Paper RSC Advances This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. Fig. 4 A centroid snapshot of the most populated cluster, chosen as a snapshot representative of the ensemble of structures generated during the last MD simulation of the Gly526Ser COX-2/AA Michaelis complex. parameters of the histidine coordinated to the metal were In the second minimization, only the heme prosthetic group described by the general AMBER force eld (GAFF)23 using the and His388 were restrained. In the last minimization only the residuegen and parmchk modules of AMBER14.24 Secondly, the protein side chains were restrained, the rest of the system AA substrate, the delocalized C11–C15 pentadienyl radical of AA, remaining free. Then, we carried out the MD simulation and its delocalized C13–C15 allyl radical parameters were also employing periodic boundary conditions (PBCs), with a cutoff taken from a recent study.25 Finally, the parameterization of the of 10 Å for all Lennard-Jones and electrostatic interactions and tyrosyl radical was performed by calculating rst the optimized the particle-mesh Ewald method31 to treat long-range electro- structure of this residue using the B3LYP functional26,27 and the static effects. The system was gradually heated from 0 to 300 K, 6-31G(d,p)28,29 basis set for all atoms, followed by an RESP with an increment of 30 K, during 10 steps of 20 ps each, calculation of the partial atomic charges at the same level of imposing a harmonic restraint (5 kcal mol 1 Å 2) to the heme theory using Gaussian09.30 At this point, as the residue is not prosthetic group, under the canonic (NVT) ensemble using a standard tyrosine, the parameters of its bonds, angles, dihe- Langevin dynamics.32 Next, a simulation was performed under drals and van der Waals terms should be re-calculated. There- the isothermic–isobaric (NPT) ensemble to ensure that the box fore, the force constants associated with each bond and angle were tted using a harmonic function to reproduce the B3LYP/ 6-31G(d,p) potential energy proles generated for the corre- sponding stretching and bending vibrations. The dihedral parameters were not calculated as they do not present signi- cant differences with the standard tyrosine residue. The van der Waals parameters were taken from the GAFF force eld of AMBER14, assigning the oxygen radical as an ether or ester oxygen. 2.3 Molecular dynamics simulations Starting from the structure resulting from the setup, we carried out three successive molecular dynamics simulations using the same protocol before the production stage. For the rst MD simulation the monomer B of Gly526Ser COX-2, with AA and the radical Tyr385 within the cyclooxygenase active site, was opti- mized following three successive molecular mechanics (MM) minimizations using the steepest descent and conjugate Fig. 5 Forward (blue line) and backward (ochre line) potential energy gradient methods. This was done to relax the system. Firstly, we profiles corresponding to C13 pro-S hydrogen abstraction by the applied harmonic restraints on the enzyme, the heme pros- tyrosyl radical. The circles indicate the locations of the stationary thetic group and AA, keeping the solvation water molecules free. points. This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 986–997 | 989 View Article Online RSC Advances Paper frozen. This region was dened by all residues and water molecules within a 15 Å radius sphere centering on C10, including the complete AA and the heme prosthetic group. The QM/MM calculations were carried out with the modular This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. package ChemShell,37,38 using Gaussian09 for the density functional theory (DFT) calculations of the quantum mechanics (QM) part. The calculations of the MM part were carried out with the DL_POLY module39 in Chemshell, employing the AMBER force eld.18 The QM region (being the Fe atom in the MM part because it did not belong to the cyclooxygenase active Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. site) was described by the B3LYP functional with the 6-31G(d,p) basis set for all atoms, with a zero total charge and doublet multiplicity in all reaction steps. This QM region presented different sizes depending on the particular step studied in the catalytic mechanism (see Fig. S1†). The polarizing effect of the Fig. 6 Potential energy profiles corresponding to the addition of the enzyme over the QM part was included by the electrostatic oxygen molecule to C11. The first point of reaction path I sets the origin embedding scheme,40 and the QM/MM boundary was treated of potential energies. O1 stands for the oxygen atom of the oxygen using hydrogen link atoms with the charge-shi model.41 molecule that is closest to C11. Finally, no cutoffs were introduced for the nonbonding MM and QM/MM interactions. To determine the potential energy proles, a series of opti- did not deform, and to achieve an adequate density. Thus, once mizations were carried out imposing harmonic restraints on the at 300 K, we carried out four steps of 10 ps each, at 1 atm reaction coordinate, with an increment of 0.2 Å at each step. pressure, with a weak restraint on the residue side chains, using The energy minimizations were done with the limited-memory an isotropic weak-coupling algorithm and the Berendsen Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm.42 A barostat,33 and then one more step of 160 ps with the same combination of the partitioned rational function optimizer (P- restraint up to a density of around 1 g cm 3. Finally, an equil- RFO)39 and the L-BFGS methods was used to locate transition ibration stage of 10 ns followed by a production period of 100 ns state structures. All these algorithms are included in the Hybrid were calculated at 300 K under the NVT ensemble without any Delocalized Internal Coordinate Optimiser (HDLCopt)39 module restraints apart from employing the SHAKE algorithm34 to of Chemshell. All the pictures of molecules were generated with constrain all the covalent bonds containing hydrogen. the VMD43 and CHIMERA17 programs. The other two MD simulations will be explained below. All the MD simulations were performed using the AMBER18 GPU (CUDA) version of the PMEMD package.35,36 3 Results and discussion 3.1 Molecular dynamics simulations 2.4 QM/MM calculations As explained above, the rst MD simulation of the Gly526Ser We used quantum mechanics/molecular mechanics (QM/MM) COX-2/AA Michaelis complex was carried out over 100 ns calculations for the energy optimization of the structures. starting from the structure resulting from the setup. The Each QM/MM system consisted of nearly 10 500 atoms, evolution of the distances between the C13 pro-S hydrogen or including all the residues of the monomer B of Gly526Ser COX-2 the C13 carbon atom and the hydrogen-acceptor oxygen of the (with the radical Tyr385), the complete AA, and the 500 water Tyr385 radical is displayed in Fig. 1 (for greater clarity, see the molecules closest to C10 of AA. The active region, where all separated evolution of both distances in Fig. S2 and S3†). The atoms moved freely during the optimization processes, involved distributions of both distances practically match throughout around 2100 atoms whereas the rest of the system was kept the simulation. The distances oscillate between about 5 Å and 7 Table 1 Initial O–C11 distances, potential energy barriers, type of oxygen addition, and stereochemistry at C11 for the six reaction paths cor- responding to the addition of the oxygen molecule to C11 Reaction path Initial distance O–C11 (Å) DE (kcal mol 1) O2 attack Stereochemistry at C11 I 3.18 7.7 Suprafacial R II 3.43 8.4 Suprafacial R III 3.12 8.6 Suprafacial R IV 4.04 11.2 Suprafacial R V 3.70 7.3 Antarafacial S VI 3.84 7.1 Antarafacial S 990 | RSC Adv., 2020, 10, 986–997 This journal is © The Royal Society of Chemistry 2020 View Article Online Paper RSC Advances Å most of the time for the rst 50 ns, but, aer a sudden jump, Table 2 Potential energy barriers, and stereochemistry at C11 and C9 both remain between about 7 Å and 9 Å for the last 50 ns. for the six reaction paths corresponding to the 9,11-cyclization leading to a C8-radical cyclic endoperoxide A typical structure corresponding to the second half of the MD simulation is pictured in Fig. 2. It can be seen that AA is Reaction DE This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. forced to undergo a major repositioning as a consequence of the path (kcal mol 1) Stereochemistry at C11 Stereochemistry at C9 Gly526Ser mutation, moving from the “L-shaped” binding conformation that it adopts in the wild-type COX-2 3 to I 8.0 R S II 7.3 R S a conformation where the chains corresponding to the AA u- III 7.7 R S end and the AA carboxylate-end become roughly parallel. This IV 8.1 R S motion causes the C13 pro-S hydrogen to be too far from the V 9.9 S R Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. hydrogen-acceptor oxygen of the Tyr385 radical for the VI 10.7 S R abstraction of the rst step of the all-radical mechanism to occur. This rather unexpected behavior of AA inside the Gly526Ser mutant could not explain the experimental results by Brash and 2.1, and we repeated the protocol described in Section 2.3 to run co-workers.16 Thus, we checked if the problem lay in the fact another 100 ns MD simulation of the Gly526Ser COX-2/AA that we generated the setup of the mutant keeping the substrate Michaelis complex. The evolution of the distances between the AA within it. That is, we had not thus far allowed a free relax- C13 pro-S hydrogen or the C13 carbon atom and the hydrogen- ation of the mutant enzyme in the absence of AA. Then, we acceptor oxygen of the Tyr385 radical (see Fig. 3) shows that eliminated AA from the setup (see Section 2.1), we repeated the most of the structures are now ready for the abstraction in the protocol previous to the production stage as described in rst step of the all-radical mechanism, especially for the rst 90 Section 2.3, and we ran a 20 ns MD simulation of the Gly526Ser ns, when the distances C13 pro-S hydrogen–O(Tyr385) and C13– COX-2 mutated enzyme. During these 20 ns, Ser526 rotates in O(Tyr385) are 3.4  0.8 Å and 3.8  0.8 Å, respectively. order that the hydrogen of its OH group forms a hydrogen bond In order to choose a snapshot representative of the ensemble with the oxygen atom of the backbone of Met522 (this hydrogen of structures generated along this last MD simulation of the bond is not possible with Gly526 in wild-type COX-2). This Gly526Ser COX-2/AA Michaelis complex, we picked one struc- position is mostly maintained during this MD simulation (see ture each 10 ps and ltered the structures by imposing the Fig. S4†). As a consequence, the two hydrogen atoms of the Cb of conditions r(C13 pro-S H–O(Tyr385)) # 3 Å and r(C13 pro-S H– Ser530 point towards the hydrophobic groove where the AA u- O(Tyr385)) < r(C13–O(Tyr385)). So, we obtained 3049 snapshots, end will nally lie. which we clustered using a root mean square deviation (RMSD) Next, we re-introduced AA into the cavity of the mutant. To of 1.3 Å for the heavy atoms of AA. The three most populated this aim, we eliminated the water molecules from the setup (see clusters included 1737, 1000 and 305 structures, respectively. Section 2.1) and from the last structure of the previous 20 ns MD Next, we determined the average position of the heavy atoms of simulation, and we overlaid the enzyme residues of both struc- AA over all the structures belonging to the most populated tures. Then, we re-solvated the system, as explained in Section cluster, and we found the snapshot of that cluster with the smallest AA heavy atoms RMSD with respect to that average position. This snapshot, the so-called centroid snapshot for the Fig. 7 Potential energy profiles corresponding to the 9,11-cyclization leading to a C8-radical cyclic endoperoxide. O2 stands for the oxygen atom of the oxygen molecule that is free and can attack C9. The Fig. 8 Potential energy profiles as a function of the C8–C12 distance potential energy of the initial structure of reaction path I is taken as the taken as the reaction coordinate. The potential energy of the initial origin of energies. structure of reaction path I is taken as the origin of energies. This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 986–997 | 991 View Article Online RSC Advances Paper most populated cluster, is pictured in Fig. 4. This structure has Fig. S5.† In addition, the distances between side chain carbon no similarities with the one shown in Fig. 2 in the case of the atoms of selected residues of that snapshot 1 (COX-2/AA) or that unrelaxed mutant enzyme. Now the “L-shaped” productive centroid snapshot (Gly526Ser COX-2/AA) and selected carbon binding conguration of AA, as observed in the above- atoms of AA bound in the corresponding active site are shown in This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. mentioned COX-2 crystal structure,3 has been recovered, with Table S1.† the AA u-end extending along the hydrophobic groove above The centroids of the second and third most populated Ser530, and the C13 carbon atom and the C13 pro-S hydrogen clusters were also determined. They are compared with the atom of AA located below and near the phenolic oxygen of centroid of the most populated cluster in Fig. S6.† An “L-sha- Tyr385. A comparison between the structure corresponding to ped” binding conguration of AA is obtained in the three cases. snapshot 1 extracted from the MD simulation of the COX-2/AA In order to check if the 100 ns MD simulation shown in Fig. 3 Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. Michaelis complex (see ref. 13) and the centroid snapshot of is actually representative of the congurational space visited by the most populated cluster of the last MD simulation of the an equilibrated Gly526Ser COX-2/AA Michaelis complex, we Gly526Ser COX-2/AA Michaelis complex (Fig. 4) is pictured in carried out two additional MD simulations. In the rst, we Fig. 9 Cyclic endoperoxide corresponding to reaction path II as a function of the C8–C12 distance at a point before the hydrogen transfer and the O–O cleavage (A), and at a point after that (B). Distances are given in Å. 992 | RSC Adv., 2020, 10, 986–997 This journal is © The Royal Society of Chemistry 2020 View Article Online Paper RSC Advances lengthened this MD simulation for an additional 150 ns (see coordinate. Only six different reaction paths leading to the Fig. S7†). In the second, we selected at random one of the peroxyl radical in C11 were obtained. In Fig. 6 we display these structures of the 20 ns MD simulation of the Gly526Ser COX-2 reactions paths (labeled I to VI) according to the structure (that mutated enzyme, and aer re-introducing AA into the cavity is, the initial optimized position of the oxygen molecule) from This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. of the mutant as explained above, we generated a 100 ns MD which the reaction path starts. The corresponding potential simulation of the Gly526Ser COX-2/AA Michaelis complex (see energy barriers and geometrical information are given in Table Fig. S8†). This new MD simulation is equivalent to the one 1. The suprafacial and antarafacial oxygen additions lead to shown in Fig. 3, but starting from a different structure and peroxides with R and S congurations at C11, respectively. This different initial velocities. Comparison of Fig. 3 and S7, S8† scenario is the opposite to the case of wild-type COX-2, where conrms that, except for some local deviations, the Gly526Ser the suprafacial and antarafacial oxygen attacks produce the S Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. COX-2/AA Michaelis complex populates a region of the cong- and R congurations at C11, respectively.13 On the other hand, urational space where most of the structures are ready for the the potential energy barriers (7.1 to 11.2 kcal mol 1) turn out to C13 pro-S hydrogen abstraction by O(Tyr385) in the rst step of be small in comparison with that for hydrogen abstraction the all-radical mechanism. (24.8 kcal mol 1). The third step consists of the 9,11-cyclization, leading to a C8-radical cyclic endoperoxide. In Fig. 7 we have plotted the 3.2 QM/MM calculations potential energy proles obtained as a function of the reaction Now we used QM/MM calculations to determine if the coordinate O–C9 for the six reaction paths. Some of the proles Gly526Ser mutant of COX-2 is capable of driving the centroid snapshot of AA throughout the all-radical mechanism. Firstly, we built the potential energy prole of the C13 pro-S hydrogen abstraction by the tyrosyl radical (step 1) as a function of the reaction coordinate. In this case, the reaction coordinate was taken as the difference between the distances corresponding to the breaking bond (C13-pro-S hydrogen) and the forming bond (pro-S hydrogen–O(Tyr385)). The energy prole was calculated both forward and backward to avoid hysteresis (see Fig. 5). Once the reactants, transition state structure, and products had been located, the potential energy barrier turned out to be 24.8 kcal mol 1. The values of the distances corresponding to the three atoms that directly participate in the breaking/ forming bonds are given in Table S2† for each structure. As expected for a hydrogen transfer, the hydrogen donor and acceptor atoms approach each other to make the hydrogen transfer easier. Aer the transfer, both atoms move away. The transition state structure appears in the region with the closest distance (see Fig. S9†). The rst step produces a planar delocalized C11–C15 penta- dienyl radical AA. In the second step the addition of the oxygen molecule to C11 can be antarafacial (with the face of the C11–C15 pentadienyl system opposite to the tyrosyl radical), giving an 11R stereochemistry, or suprafacial (on the same side as the tyrosyl radical), leading to an 11S stereochemistry. We chose different possible starting positions of the oxygen molecule for the attack to C11, as we had already done in previous work.13 That C11 atom was chosen as the origin of coordinates, and the oxygen molecules, when possible, were placed around it along the x, y and z Cartesian axes, along the bisector axes contained in the xy, xz and yz planes, or near them. In all, 53 initial oxygen molecules were set close to C11, at distances of 2.5, 3, and 3.5 Å, corresponding to 9, 18, and 26 O2 molecules, respectively. Fiy- three QM/MM single point energy calculations were carried out for each one. Then, the structures with the highest energies Fig. 10 Potential energy profiles corresponding to the formation of (A) the 8,9- and 11,12-epoxides in a concerted way for reaction paths I and were discarded and the most stable structures were optimized. II, but just the 8,9-epoxide for reaction paths III–VI, as a function of the The optimized structures were taken as starting points to build O2 (bonded to C9)–C8 distance taken as the reaction coordinate; and the potential energy proles for the oxygen addition to C11, the (B) the 11,12-epoxide for reaction paths III to VI, as a function of the O1 distance O (the one closest to C11)–C11 being the reaction (bonded to C11)–C12 distance taken as the reaction coordinate. This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 986–997 | 993 View Article Online RSC Advances Paper involve two potential energy maxima. In these cases, in the rst 11R) and (9R, 11S) cyclic endoperoxides, but not (9R, 11R). This stage the free oxygen of the oxygen molecule added to C11 has to fact will have important consequences for the following mech- rotate around the C11–O bond in order to approach C9 (rst anistic step. maximum) and then the 9,11-cyclization occurs (second For the fourth step we chose the distance C8–C12 to dene This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. maximum). The potential energy barriers corresponding to the the reaction coordinate for the 8,12-cyclization. When that higher energy transition state structures for each reaction path distance was shortened the O–O bond of the cyclic endoper- are shown in Table 2, along with stereochemical information. oxide was broken in all cases. Instead of the formation of the The four paths corresponding to an R conguration at C11 cyclopentane ring, a hydrogen transfer from C11 to C8 followed involve barriers 2–3 kcal mol 1 lower than the two paths with an by the breakage of the O–O bond occurred, leading to the S conguration at C11. Most interestingly, we have obtained (9S, formation of a ketone group in C11 and a radical centered at the Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. Fig. 11 Structure of the 8,9- and 11,12-epoxides corresponding to (A) reaction path I; and (B) reaction path III. Distances are given in Å. 994 | RSC Adv., 2020, 10, 986–997 This journal is © The Royal Society of Chemistry 2020 View Article Online Paper RSC Advances other oxygen atom (the one bonded to C9). The corresponding epoxide in the case of reaction paths III and IV (Fig. 10B) just potential energy proles are displayed in Fig. 8. This process correspond to the previous reorganization of Tyr385 in order to turns out to be very exothermic, and the hydrogen transfer and be able to interact with the epoxide to be formed. It is inter- the O–O cleavage can take place in a concerted way or in two esting to note that, conversely, Tyr385 does not interact with This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. stages, depending on the reaction path. The potential energy either of the two oxygen atoms of the O–O bridge aer the 8,12- barriers range from 15.5 kcal mol 1 to 25 kcal mol 1. cyclization to give a bicyclo endoperoxide which occurs in wild- Our calculations conrm the experimental nding by Brash type COX-2 (see Fig. S10,† derived from the calculations in our and co-workers,16 who found that the 8,12-cyclization does not previous work13 for the sake of comparison). occur in the Gly526Ser mutant COX-2. We know that the cyclopentane ring closure only occurs with a (9R, 11R) cong- Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. uration,13,16 but this is just the stereochemical conguration 4. Conclusions that cannot be reached aer the third mechanistic step in the mutant, as mentioned above. The reason for that can be The biosynthesis of prostaglandins from arachidonic acid understood by analyzing Fig. 9, where, for the sake of example, catalyzed by cyclooxygenase-2 is the result of a complex mech- we have pictured the cyclic endoperoxide corresponding to anism where every residue of the cyclooxygenase active site reaction path II as a function of the C8–C12 distance at a point plays its role to achieve exceptional stereochemical and regio- before the hydrogen transfer and the O–O cleavage (Fig. 9A), chemical control of the outcome. Unraveling these roles is not and at a point aer that (Fig. 9B). As commented in Section 3.1, an easy task. Particularly remarkable is the case of the Gly526 in the Gly526Ser mutant COX-2 the hydrogen of the Ser526 OH residue. Brash and co-workers16 have shown that the Gly526Ser group forms a hydrogen bond (see Fig. S5†) with the oxygen mutant of human COX-2 is unable to produce the 8,12-cycliza- atom of the Met522 backbone (a hydrogen bond that cannot tion required for the formation of prostaglandin. Instead, they exist with Gly526 in wild-type COX-2). As a consequence, the two found novel products that were identied as 8,9-11,12-diepoxy- hydrogen atoms of the Cb Ser530 point towards the hydrophobic 13R- (or 15R)-hydroperoxy derivatives of arachidonic acid. In groove where the AA u-end lies. This position makes the helix D this work we have combined molecular dynamics simulations containing residues from 521 to 535 quite rigid and forces and QM/MM calculations to shed light on the molecular details Ser530 to adopt the position shown in Fig. 9A (both Ser residues of the effects of such a mutation. belong to that helix). As a result, the hydrogen at C9 cannot be at Our results show that inside the Gly526Ser COX-2 binding the same side as the hydrogen at C11 (this is equivalent to saying pocket, AA adopts the same “L-shaped” productive binding that the (9R, 11R) conguration cannot be reached), and the conguration as the one observed in the COX-2 crystal struc- 8,12-cyclization cannot happen. In Fig. 9B the hydrogen transfer ture,3 with the AA u-end extending along the hydrophobic from C11 to C8 and the breakage of the O–O bond have already groove above Ser530, and the C13 carbon atom and the C13 pro-S taken place. hydrogen atom of AA located below and near the phenolic At this point it is clear that the Gly526Ser mutation hinders oxygen of Tyr385. Then, the rst three reaction steps of the all- the formation of prostaglandins because the 8,12-cyclization radical mechanism7,13 can take place, leading to C8-radical does not occur, and the cyclic endoperoxide breaks. Neverthe- cyclic endoperoxides, in both wild-type COX-2 and Gly526Ser less, this breakage does not experimentally take place in the way COX-2. The difference between the two cases lies in the fact that described above according to Fig. 8 and 9B. In effect, as already only COX-2, but not the mutant, is able to produce cyclic found by Brash and co-workers,16 the homolytic cleavage of the endoperoxides with a (9R, 11R) conguration, the stereochem- O–O bond actually leads to the formation of 8,9-11,12-diepoxy istry that is required to make the 8,12-cyclization possible. This derivatives of the arachidonic acid. Taking the distance O is a consequence of the fact that in the Gly526Ser mutant COX-2, (bonded to C9)–C8 as the reaction coordinate, the 8,9- and 11,12- the hydrogen of the Ser526 OH group forms a hydrogen bond epoxides are formed in a concerted way for reaction paths I and with the oxygen atom of the Met522 backbone (a hydrogen bond II (see Fig. 10A), but just the 8,9-epoxide for reaction paths III–VI that cannot exist with Gly526 in wild-type COX-2). Then, the two (see also Fig. 10A). The 11,12-epoxide for these last reaction hydrogen atoms of Cb Ser530 point towards the hydrophobic paths is formed next when the distance O (bonded to C11)–C12 is groove where the AA u-end is located, in a position that hinders chosen to dene the reaction coordinate (see Fig. 10B). Anyway, the hydrogen at C9 from being at the same side as the hydrogen the product contains an allyl radical delocalized from C13 to C15. at C11 (this is equivalent to saying that the (9R, 11R) congu- In all cases the potential energy barriers for the formation of the ration cannot be reached). Consequently, the 8,12-cyclization epoxides are clearly lower than the ones corresponding to the towards the formation of prostaglandin cannot take place in above-discussed O–O cleavage involving a hydrogen transfer Gly526Ser COX-2. In fact, all attempts to force the shortening of from C11 to C8, so explaining the experimental formation of the the C8–C12 distance provoke hydrogen transfer from C11 to C8 8,9-11,12-diepoxy derivatives. followed by the breakage of the O–O bond of the cyclic endo- Both epoxides are stabilized by interactions with several peroxide. Anyway, the homolytic cleavage of the O–O bond residues like Ser526. The 11,12-epoxide is especially stabilized leading to the formation of 8,9-11,12-diepoxy derivatives of the by formation of a hydrogen bond with Tyr385, which, in turn, is arachidonic acid turns out to occur faster, thus explaining the hydrogen-bonded to Tyr348 (see Fig. 11). As a matter of fact, the experimental ndings by Brash and co-workers.16 We hope that two early maxima that appear for the formation of the 11,12- this work can contribute to better molecular comprehension of This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 986–997 | 995 View Article Online RSC Advances Paper the catalytic mechanism of COX-2, one of the main enzymes A. W. Goetz, D. Greene, R. Harris, N. Homeyer, S. Izadi, responsible for inammation processes in humans. A. Kovalenko, T. Kurtzman, T. S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, D. J. Mermelstein, Conflicts of interest K. M. Merz, Y. Miao, G. Monard, C. Nguyen, H. Nguyen, This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. I. Omelyan, A. Onufriev, F. Pan, R. Qi, D. R. Roe, There are no conicts to declare. A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C. L. Simmerling, J. Smith, R. Salomon-Ferrer, J. Swails, Acknowledgements R. C. Walker, J. Wang, H. Wei, R. M. Wolf, X. Wu, L. Xiao, D. M. York and P. A. Kollman, AMBER 2018, University of This work was supported by the Spanish “Ministerio de Ciencia, California, San Francisco, 2018. Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. Innovación y Universidades” (project CTQ2017-83745-P). 20 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926– References 935. 21 K. Shahrokh, A. Orendt, G. S. Yost and T. E. Cheatham, J. 1 M. Murakami, Exp. Anim., 2011, 60, 7–20. Comput. Chem., 2012, 33, 119–133. 2 C. Schneider, D. A. Pratt, N. A. Porter and A. R. Brash, Chem. 22 C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Biol., 2007, 14, 473–488. Chem., 1993, 97, 10269–10280. 3 A. J. Vecchio, D. M. Simmons and M. G. Malkowski, J. Biol. 23 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and Chem., 2010, 285, 22152–22163. D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174. 4 C. Yuan, C. J. Rieke, G. Rimon, B. A. Wingerd and 24 D. A. Case, V. Babin, J. T. Berryman, R. M. Betz, Q. Cai, W. L. Smith, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 6142– D. S. Cerutti, T. E. Cheatham III, T. A. Darden, R. E. Duke, 6147. H. Gohlke, A. W. Goetz, S. Gusarov, N. Homeyer, 5 W. L. Smith, D. L. DeWitt and R. M. Garavito, Annu. Rev. P. Janowski, J. Kaus, I. Kolossváry, A. Kovalenko, T. S. Lee, Biochem., 2000, 69, 145–182. S. LeGrand, T. Luchko, R. Luo, B. Madej, K. M. Merz, 6 B. J. Orlando, P. P. Borbat, E. R. Georgieva, J. H. Freed and F. Paesani, D. R. Roe, A. Roitberg, C. Sagui, R. Salomon- M. G. Malkowski, Biochemistry, 2015, 54, 7309–7312. Ferrer, G. Seabra, C. L. Simmerling, W. Smith, J. Swails, 7 M. Hamberg and B. Samuelsson, J. Biol. Chem., 1967, 242, R. C. Walker, J. Wang, R. M. Wolf, X. Wu and 5336–5343. P. A. Kollman, AMBER 14, University of California, San 8 K. Marnett and L. J. Maddipati, ed. M. Grisham, J. Everse and Francisco, 2014. K. Everse, Prostaglandin H Synthase, CRC Press, Florida, 25 P. Tosco, J. Am. Chem. Soc., 2013, 135, 10404–10410. Unites States, Boca Raton, 1991, vol. 1, pp. 293–334. 26 C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter 9 W. L. Smith and L. J. Marnett, Biochim. Biophys. Acta, 1991, Mater. Phys., 1988, 37, 785–789. 1083, 1–17. 27 A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652. 10 C. Schneider, W. E. Boeglin, S. Lai, J. K. Cha and A. R. Brash, 28 P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, Anal. Biochem., 2000, 284, 125–135. 213–222. 11 L. J. Marnett, Curr. Opin. Chem. Biol., 2000, 4, 545–552. 29 M. M. Francl, W. J. Pietro, W. J. Hehre, J. S. Binkley, 12 P. J. Silva, P. A. Fernandes and M. J. Ramos, Theor. Chem. M. S. Gordon, D. J. DeFrees and J. A. Pople, J. Chem. Phys., Acc., 2003, 110, 345–351. 1982, 77, 3654–3665. 13 A. Cebrián-Prats, À. González-Lafont and J. M. Lluch, ACS 30 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, Omega, 2019, 4, 2063–2074. M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, 14 S. W. Rowlinson, B. C. Crews, D. C. Goodwin, C. Schneider, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, J. K. Gierse and L. J. Marnett, J. Biol. Chem., 2000, 275, 6586– X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, 6591. J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, 15 C. Schneider, W. E. Boeglin, J. J. Prusakiewicz, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, S. W. Rowlinson, L. J. Marnett, N. Samel and A. R. Brash, J. O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Biol. Chem., 2002, 277, 478–485. J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, 16 C. Schneider, W. E. Boeglin and A. R. Brash, J. Biol. Chem., E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, 2004, 279, 4404–4414. J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, 17 E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, D. M. Greenblatt, E. C. Meng and T. E. Ferrin, J. Comput. M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, Chem., 2004, 25, 1605–1612. J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, 18 J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, K. E. Hauser and C. Simmerling, J. Chem. Theory Comput., R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, 2015, 11, 3696–3713. P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, 19 D. A. Case, I. Y. Ben-Shalom, S. R. Brozell, D. S. Cerutti, Farkas, J. B. Foresman, J. V Ortiz, J. Cioslowski, D. J. Fox, T. E. Cheatham III, V. W. D. Cruzeiro, T. A. Darden, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, R. E. Duke, D. Ghoreishi, M. K. Gilson, H. Gohlke, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, 996 | RSC Adv., 2020, 10, 986–997 This journal is © The Royal Society of Chemistry 2020 View Article Online Paper RSC Advances B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, 37 S. Metz, J. Kästner, A. A. Sokol, T. W. Keal and P. Sherwood, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 101–110. J. L. Sonnenberg, M. Had and D. J. Fox, Gaussian, Inc., 38 P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, Wallingford CT, 2009. C. R. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. 31 M. Deserno and C. Holm, J. Chem. Phys., 1998, 109, 7678– W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, 7693. J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, 32 H. Lei, N. A. Baker and X. Li, Proc. Natl. Acad. Sci. U. S. A., E. Karlsen, M. Sjøvoll, A. Fahmi, A. Schäfer and 2016, 113, 14183–14188. C. Lennartz, J. Mol. Struct.: THEOCHEM, 2003, 632, 1–28. 33 H. J. C. Berendsen, J. P. M. Postma, W. F. Van Gunsteren, 39 W. Smith and T. R. R. Forester, J. Mol. Graphics, 1996, 14, A. Dinola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690. 136–141. Open Access Article. Published on 03 January 2020. Downloaded on 6/9/2020 11:31:35 AM. 34 J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. 40 D. Bakowies and W. Thiel, J. Phys. Chem., 1996, 100, 10580– Phys., 1977, 23, 327–341. 10594. 35 S. Le Grand, A. W. Götz and R. C. Walker, Comput. Phys. 41 H. M. Senn and W. Thiel, Angew. Chem., Int. Ed., 2009, 48, Commun., 2013, 184, 374–380. 1198–1229. 36 R. Salomon-Ferrer, A. W. Götz, D. Poole, S. Le Grand and 42 D. C. Liu and J. Nocedal, Math. Program., 1989, 45, 503–528. R. C. Walker, J. Chem. Theory Comput., 2013, 9, 3878–3888. 43 W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38. This journal is © The Royal Society of Chemistry 2020 RSC Adv., 2020, 10, 986–997 | 997