Molecular mechanisms in the selectivity of non-steroidal anti

Jan 18, 2018 - Nonsteroidal anti-inflammatory drugs (NSAIDs) inhibit cyclooxygenase (COX) 1 and 2 with varying degrees of selectivity. A group of COX-...
0 downloads 0 Views 7MB Size
Article Cite This: Biochemistry 2018, 57, 1236−1248

pubs.acs.org/biochemistry

Molecular Mechanisms in the Selectivity of Nonsteroidal AntiInflammatory Drugs Yasmin Shamsudin Khan, Hugo Gutiérrez-de-Terán, and Johan Åqvist* Department of Cell and Molecular Biology, Box 596, Uppsala University, BMC, SE-751 24 Uppsala, Sweden S Supporting Information *

ABSTRACT: Nonsteroidal anti-inflammatory drugs (NSAIDs) inhibit cyclooxygenase (COX) 1 and 2 with varying degrees of selectivity. A group of COX-2 selective inhibitorscoxibsbinds in a time-dependent manner through a three-step mechanism, utilizing a side pocket in the binding site. Coxibs have been extensively probed to identify the structural features regulating the slow tightbinding mechanism responsible for COX-2 selectivity. In this study, we further probe a structurally and kinetically diverse data set of COX inhibitors in COX-2 by molecular dynamics and free energy simulations. We find that the features regulating the high affinities associated with time-dependency in COX depend on the inhibitor kinetics. In particular, most time-dependent inhibitors share a common structural binding mechanism, involving an induced-fit rotation of the side-chain of Leu531 in the main binding pocket. The high affinities of two-step slow tight-binding inhibitors and some slow reversible inhibitors can thus be explained by the increased space in the main binding pocket after this rotation. Coxibs that belong to a separate class of slow tight-binding inhibitors benefit more from the displacement of the neighboring side-chain of Arg513, exclusive to the COX-2 side-pocket. This displacement further stabilizes the aforementioned rotation of Leu531 and can explain the selectivity of coxibs for COX-2.

P

COX-1, which is Val523 in COX-2. Since this residue is positioned between the main and side pockets (Figure 1), it was originally believed that the bulkier Ile523 in COX-1 would block access to the side-pocket. The selectivity of coxibs toward COX-2 would thus be caused by the increased facility to access the side-pocket of COX-2.13 However, it has since been demonstrated that coxibs can bind in the corresponding pocket in COX-1.14 An alternative explanation for selectivity is based on the kinetic profile of COX inhibitors.15 NSAIDs can be divided into five kinetic classes, which differ in the reversibility and number of steps involved in the binding mechanism.16 Rapidly reversible inhibitors, such as ibuprofen (Figure 2), bind according to the simple one-step mechanism

ain-relieving and fever-reducing nonsteroidal anti-inflammatory drugs (NSAIDs) inhibit both isoforms of cyclooxygenase (COX). COX-1 is constitutively expressed,1 while COX-2 is induced in response to inflammation.2 Most traditional NSAIDs inhibit both isoforms in a nonselective way, but it is believed that it is the inhibition of COX-2 that results in pain relief.3 Moreover, equally or more potent inhibition of COX-1 is responsible for the common side-effects of NSAIDs, such as peptic ulcers, gastrointestinal bleeding,4 and kidney problems.5,6 Consequently, there is a strong focus on the development of pain-killer drugs that selectively target COX-2, which resulted in the release of a new class of NSAIDs, often referred to as “coxibs”. However, some of these were withdrawn from the market due to their increased risk of myocardiac infarctions.7 Still, coxibs are commonly used to treat arthritis, chronic pain,8 and colon cancers.9 In the cancer area, the selective inhibition of COX-2 is a promising drug target, since it was found that this enzyme isoform is overexpressed in many different human tumors cells,10,11 but not the surrounding healthy cells.12 However, the underlying structural reasons for COX selectivity have remained elusive. The two isoforms share the same tertiary structure with 58% sequence identity. To facilitate comparison between the two isoforms the labeling of residue numbers in hCOX-2 in this work will follow that of the crystal structures of COX-1. Several points of variability are located within the binding pocket and the adjacent side-pocket (Figure 1), with some having been explicitly singled out as selectivity hotspots. One of these is Ile523 in the main binding pocket of © 2018 American Chemical Society

E + I ⇄ EI

(1)

The only known truly irreversible NSAID is acetylsalicylic acid (aspirin) (Figure 2). It binds irreversibly in a time-dependent manner

E + I ⇄ EI → EI*

(2)

In the first step the ligand binds to the reactive site, and this binding event is reversible. In the second step Ser530 in the binding pocket is acetylated, which results in a noncovalently bound side-product of salicylic acid.17 Received: October 9, 2017 Revised: January 8, 2018 Published: January 18, 2018 1236

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

Figure 1. Detailed view of the binding sites and side-pockets of the generated homology models of human COX-1 (gray) and COX-2 (magenta). Identical residues in both isoforms are colored green. Leu531 is shown in both the open (blue) and closed (yellow) rotational states.

slow tight-binding inhibitors, but binds through a three-step mechanism E + I ⇄ EI ⇄ EI* → EX

where the third step is time-dependent and only slowly reversible. These inhibitors, which include the coxib drug celecoxib (Figure 2), are time-dependent in COX-2 but not in COX-1.15 Because of this different kinetic profile for each COX isoform, it has been suggested that the time-dependence causes selectivity. However, the structural features of this three-step mechanism have yet to be identified, which is the focus of the present study. Crystal structures of ovine COX-1 and murine COX-2 show that COX has at least two stable conformations, mainly identified by the rotational conformation of the Leu531 sidechain.18,19 The most common conformation observed in crystal structures, referred to as the “closed” conformation, presents the Leu531 side-chain pointing into the binding pocket (Figure 1). Through conformational analyses of COX-1 we have previously demonstrated that this is the preferred conformation in apo human COX-1 (hCOX-1) and also in complex with the rapidly reversible inhibitor ibuprofen (Figure 2).20,21 However, two-step slow tight-binding inhibitors (eq 2) stabilize a secondary “open” conformation, which in turn induces higher binding affinities of the inhibitors. Herein, we extend our analysis to the computational characterization of ligand binding in COX-2. Our approach includes homology modeling, molecular dynamics (MD), binding free energy calculations, and potential of mean force (PMF) simulations of the COX enzymes.20,22 The molecular mechanism behind different kinetic profiles of COX inhibitors is thus elucidated, including the three-step mechanism of slow tight-binding inhibition of coxibs, which is exclusive to COX-2. Furthermore, to provide a deeper exploration of the selectivity hotspots, we calculated the binding affinity change upon the V523I mutation for different classes of inhibitors, through free energy perturbation (FEP) calculations. The structural rationalization of COX selectivity provided here should thus be useful in the design of new COX2 inhibitors.

Figure 2. Chemical structures of the inhibitor subset reported in Table 1 and relevant ones from Table 2. The names of the structures are color-coded according to kinetic class: rapidly reversible (orange), irreversible (green), slowly reversible (blue), two-step slow tightbinding (purple), three-step slow tight binding (red), and unknown (gray).

In between these two extreme cases lie the three classes of time-dependent NSAIDs. These include the slowly reversible, which binds according to

E + I ⇄ EI ⇄ EI*

(4)

(3)

where both steps are reversible. The next class is the slow tightbinding inhibitors where, despite the fact that they are technically reversible, the final step is only very slowly reversible, making them behave like irreversible inhibitors under physiological conditions. Thus, in practice slow tightbinding inhibitors bind in a two-step mechanism according to eq 2, with the only difference that the second step is timedependent and slowly reversible. Finally, the fifth class is also



METHODS Homology Modeling and Molecular Docking. Human cyclo-oxygenase structures were created using default homol1237

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry ogy modeling routines in Modeler,23 using an ovine COX-124 and a murine COX-2 structure25 as the templates for hCOX-1 and hCOX-2 (PDB codes 1Q4G and 3NT1 respectively). The two structures were superimposed onto the ovine COX-1 structure, and the heme group was imported from the template crystal structures by structural superposition. Hydrogen atoms were added and the structure was energy minimized using Macromodel with the OPLS_2005 force field. 26 The conformation of these original structures, similar to the template crystal structures, are hereafter referred to as the closed conformation. The open conformations, which are the induced secondary conformations, were created by rotating the Leu531 χ1 dihedral angles.18,20 This was done in the apo structure using the methodology described in the PMF section below. The average structure of the final window was subjected to a short energy minimization run and cleared of water molecules. The resulting structure is thus very similar to the open conformation structure from the PMF calculations, and also to the open conformation crystal structure of murine COX-2.18 All of the 37 COX inhibitors probed were built using the software Maestro included in the Schrödinger Suite 2011 package.27 Carboxylated ligands were modeled as negatively charged, and for chiral compounds only the S-enantiomers were modeled. To obtain initial poses for the molecular dynamics (MD) simulations all ligands were docked to both the open and closed conformations using GLIDE XP.28 Ligands were docked inside a box with the dimensions 20 × 10 × 20 Å centered on the Val523 residue. One to three poses were generated per inhibitor, and all were retained for MD simulations. An additional pose was retrieved by superimposing the energetically optimal structures of identical ligands previously docked (automated or manually) in hCOX-1. Average structures of these ligands after docking and MD simulations have been shown to be in good agreement with available crystal structures,22 which are very similar for COX-1 and COX-2 complexes. The force field parameters needed for the MD simulations of the ligands were retrieved from automatic parametrization performed with Macromodel.26 MD Simulations. Docking was followed by solvation using TIP3P29 water and MD simulations in both open and closed conformations. All MD simulations were performed using spherical boundary conditions and the OPLS-AA force field30 with the MD program Q.31 A sphere of 30 Å radius centered on Cβ atom of Ser530 was solvated with TIP3P water molecules and subject to polarization and radial constraints according to the surface constrained all-atom solvent (SCAAS) model31,32 at the sphere surface to mimic the properties of bulk water. Protein atoms outside the simulation sphere were restrained to their initial positions and only interacted with the system through bonds, angles, and torsions. All titratable residues within 20 Å of the sphere center were treated as ionized, while the remaining residues inside the simulation sphere were manually assessed with the following residues treated as ionized: Arg433, Glu465, Lys468, Lys473, Glu480, Glu502, Lys511. Titratable residues close to the sphere boundary were modeled in their neutral form to account for dielectric screening. With this setup the simulation sphere was overall neutral, except for the charge of ligands bearing a carboxylate group. Thus, the protein and water (reference) simulation systems have the same net charge, and the consideration of additional Born terms31 in the calculation of free energies is avoided.

All MD simulations started with a heating and equilibration phase followed by subsequent data collection. The system was thus gradually heated to a target temperature of 298 K, and positional restraints on all solute heavy atoms were gradually released. An MD time step of 2 fs was used for the data collection phase, and water bonds and angles, as well as solute bonds, were constrained using the SHAKE33 algorithm. In the water (reference state) simulations a weak harmonic restraint was applied to the center of mass of the ligands to keep them centered in the water droplet. Nonbonded interactions were calculated explicitly up to a 10 Å cutoff, except for the ligand atoms, for which no cutoff was used. Beyond the cutoff, longrange electrostatics were treated with the local reaction field multipole expansion method.34 Nonbonded pair lists were updated every 25 steps, and the same interval was used for the sampling of the ligand-surrounding interaction energies. Each ligand was subjected to three MD simulations in the reference state where each heating and equilibration scheme lasted 2 ns, followed by 5 ns of data collection. In an initial screening every ligand−enzyme complex pose was subjected to one MD simulation using the same heating and equilibration scheme as in the reference state simulations, followed by 10 ns of data collection. Each pose was then scored using the linear interaction energy (LIE) method. The initial pose of each ligand resulting in the lowest relative LIE binding free energy from this screening was selected as the final starting pose and was subsequently subjected to five independent MD simulations using a scheme consisting of 3 ns of heating, 2 ns unrestrained equilibration, and 10 ns of data collection. These simulation replicas were initiated with different random velocities, but otherwise set up with identical conditions. The convergence of the calculated binding free energies was judged both on the basis of the convergence of the energy components from each individual simulation replica and in terms of the final standard errors or the mean (s.e.m.) for all replicas. The individual simulations were thus considered converged if the average polar and nonpolar interaction energies had an error bar < ±1 kcal/mol based on splitting the trajectory in two halves. For the final calculated binding free energies the average s.e.m. over all 37 inhibitors is 0.51 kcal/mol, and for the subset of 13 inhibitors studied in more detail the corresponding value is 0.42 kcal/mol. Although there are a few low affinity ligands with an s.e.m. < 1.5 kcal/mol, the general convergence of the calculated binding free energies can thus be considered quite satisfactory. Potential of Mean Force Calculations. Potential of mean force (PMF) calculations were performed using the Cγ−Cβ− Cα−N dihedral angle (χ1) of Leu531 as a variable. Following 10 ns of MD simulations in the closed conformation according to the procedure described above and an additional 10 ps equilibration phase where the MD step size was reset to 1 fs, the dihedral angle was pushed from 90 to 190 deg in 100 discrete windows, each of 30 ps duration.20,35 Here, the transformation from the closed to the open state utilizes the biased force field potentials FF U O = U FF − Utor (χ ) + U bO(χ )

FF U C = U FF − Utor (χ ) + U bC(χ )

(5)

UFF tor(χ)

for the open and closed states, respectively. is the torsion term for χ1 in the regular force field potential (UFF) and UOb and UCb are biasing potentials for the open and closed states, respectively. These are quadratic, Ub = k(χ − χ0)2, with a 75 1238

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

Figure 3. (A) Average free energy profiles from PMF calculations of hCOX-2 in the apo form (blue) and in complex with the irreversible inhibitor acetylsalicylic acid (aspirin, green), rapidly reversible inhibitor ibuprofen (orange), and the slow tight-binding inhibitor flurbiprofen (purple). (B) Detailed view of the binding site and side-pocket of COX-2. The closed conformation (blue) has Leu531 pointing toward the binding pocket and a shorter distance between Arg513 and Arg120 compared to the open conformation (beige) in which Leu531 is pointing away from the main pocket.

kcal mol−1 rad−2 force constant (k) and have their respective minima (χ0) coinciding with the open and closed conformations. The transformation is then performed via a free energy perturbation mapping potential Um = (1 − λm)UO + λmUC, where the mapping parameter λm is varied between 0 and 1 in 100 discrete simulations windows. The free energy change on the true force field potential surface is then obtained by removing the bias according to the umbrella sampling formula ΔG(X n) =



in the two states referred to above, i.e., in water and embedded in the active site of the solvated protein. The nonpolar ⟨Uvdw l−s ⟩ and polar ⟨Uell−s⟩ interaction energies are treated separately, and their difference is scaled by the scaling factors α and β. The constant α has been empirically determined as 0.18, while β varies depending on the properties of the ligand.37,38 In this data set the value of β is 0.33 for neutral ligands with more than two hydroxyl groups, 0.37 for neutral ligands with one hydroxyl group, 0.43 for all remaining neutral ligands, and 0.5 for all ligands bearing a negatively charged carboxylate group. The offset term γ does not affect relative binding free energies and was thus treated as a free parameter and empirically adjusted to minimize the absolute error between calculated and experimental binding free energies, leading to γ = −4.9 kcal/mol. An electrostatic correction term39 ΔGelcorr was calculated for the interactions of charged ligands with distant neglected ionized groups neutralized in the simulations (in kcal/mol) qpql el ΔGcorr = 332 ∑ ∑ εrpl p l (9)

wm[ΔG(λm) −

m ⊂ Xn

RT ln⟨e−[U

FF

(X n) − Um(X n)]/ RT

⟩m ]/

∑ m ⊂ Xn

wm (6)

where Xn is our discretized reaction coordinate, defined as Xn = UO − UC, and the relevant energy values (UFF and Um) are binned with respect to Xn. The first term on the right-hand side of eq 6 is the free energy associated with moving on the mapping potential and is computed as the average of forward and reverse application of the free energy perturbation equation m−1

ΔG(λm) = −RT

∑ ln⟨e−(U

n + 1− Un)/ RT

⟩n

n=0

Here, qp is the formal charge of the neglected ionic group while ql is the partial charge of the ligand atom, and the effective dielectric constant ε was set to 80. The reported experimental binding free energies (ΔGexp bind) are average values of several assays, and were calculated from IC50 values according to

(7)

Since a given value of reaction coordinate Xn is sampled by multiple windows, their contribution is weighted in eq 6 with the number of data points from each λm by the factor wm/Σwm. The PMF simulations were repeated five times for each system, and the reported free energy profile is the average over these five replicas. To judge convergence, standard errors of the mean (s.e.m.) were calculated over all pooled data points in each bin, and the maximum s.e.m. was 0.75 kcal mol−1. The average hysteresis error obtained from forward and reverse application of eq 7 for the entire transformation is −4.5 meclofenamate −7.1 ± 0.1 Slow tight-binding inhibitors: E + I ⇄ EI ⇄ EI* → EX celecoxib

ΔGcalc(closed) bind

ΔGcalc(open) bind

−7.1 ± 0.2

−7.4 ± 0.2

−7.1 ± 0.2

−8.4 ± 0.5

−7.0 −9.0 −9.4 −6.6 −9.3 −8.6

0.4 0.6 n.d. 0.1 0.2 0.5

−7.0 −9.0 −7.5 −6.4 −9.4 −5.8

± ± ± ± ± ±

0.2 0.2 0.2 0.6 0.2 1.5

−8.7 − −9.4 −7.2 − −8.9

−11.1 ± 0.1 −8.2 ± 0.3 −9.7 ± 0.2 −9.5 ± 0.0

−6.7 −7.8 −5.3 −8.6

± ± ± ±

0.4 0.5 0.4 0.4

−10.9 ± 0.3 −8.8 ± 0.5 −9.1 ± 0.5 −9.8 ± 0.5

−9.2 ± n.d.

Table 3

± ± ± ± ± ±

± 0.4 ± 1.0 ± 0.2 ± 0.6

−9.2 ± 0.1

All energies are in kcal mol−1, with error bars ±1 s.e.m. for calculated values. γ = −4.9. bBinding free energies from experimental instant inhibition assays ±1 s.d.52 cAverage binding free energies from preincubated assays ±1 s.e.m.47−51,62 n.d., not determined or reported. Ligands for which the open conformation is unstable are denoted by ′−′. a

Table 2. Calculated and Experimental Binding Free Energies for Inhibitors in Complex with hCOX-2a ligand ampyrone DFP diflunisal etodolac fenoprofen ketoprofen ketorolac L745 nabumetone niflumic acid NS398 paracetamol rofecoxib salicylaldehyde salicylic acid SC58125 sulindac sulindac sulfide suprofen tenidap tolmetin tomoxiprol valerylsalicylate zomepirac

⟨Uvdw l‑s ⟩reference −19.3 −27.0 −12.2 −17.0 −14.7 −15.8 −15.7 −29.0 −23.8 −14.2 −28.4 −12.8 −27.0 −11.1 −4.2 −28.2 −21.9 −22.2 −16.0 −26.8 −16.1 −31.6 −13.0 −17.3

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0

⟨Uell‑s⟩reference −21.6 −36.9 −160.8 −191.7 −173.7 −175.7 −177.0 −31.6 −17.7 −164.8 −26.9 −34.1 −37.8 −18.0 −159.2 −35.4 −187.3 −173.0 −176.6 −35.4 −178.8 −26.2 −172.9 −175.9

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.0 0.0 0.1 0.3 0.1 0.2 0.1 0.0 0.0 0.1 0.2 0.0 0.1 0.0 0.1 0.2 0.2 0.1 0.2 0.1 0.1 0.0 0.0 0.2

⟨Uvdw l‑s ⟩complex −33.2 −46.3 −22.6 −29.7 −28.2 −29.9 −31.3 −45.4 −38.2 −26.7 −43.7 −23.4 −45.2 −19.6 −10.3 −47.1 −36.4 −43.0 −30.5 −41.8 −30.4 −49.7 −23.9 −32.0

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.7 0.7 0.6 2.3 2.3 0.1 0.5 0.2 1.3 0.4 0.6 0.5 0.7 0.1 0.3 1.2 1.0 1.2 0.3 1.5 0.2 0.2 0.6 0.8

⟨Uell‑s⟩complex −18.4 −33.9 −159.6 −193.2 −174.4 −177.1 −179.8 −28.4 −12.0 −163.0 −25.3 −29.2 −37.1 −16.7 −162.5 −35.4 −185.7 −172.5 −176.2 −35.0 −180.3 −24.3 −177.8 −175.3

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

1.0 0.3 1.2 0.6 1.7 2.8 0.0 0.9 0.5 4.8 1.2 0.5 0.1 0.7 1.0 0.8 1.9 1.2 1.2 4.1 0.5 0.2 0.2 0.7

b ΔGcalc bind

−6.0 −7.1 −6.8 −8.5 −8.3 −8.7 −9.8 −6.5 −5.0 −6.9 −7.0 −5.0 −7.9 −6.0 −8.3 −8.3 −7.4 −9.0 −8.0 −7.4 −8.9 −7.4 −9.9 −7.9

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.5 0.2 0.5 0.5 1.1 1.2 0.1 0.4 0.3 2.1 0.5 0.2 0.1 0.3 0.4 0.4 0.9 0.6 0.5 1.8 0.2 0.1 0.0 0.4

lowest energy conformation closed open closed closed open closed open open closed open open closed open closed closed open open open closed closed closed open open open

c ΔGexp(inc.) bind

−5.8 −9.6 −5.5 −8.3 −7.2 −9.1 −9.8 −8.2 −5.8 −7.0 −9.7 −6.5 −9.2 −4.7 −6.9 −9.0 −6.3 −8.2 −7.2 −7.8 −8.2 −9.2 −10.5 −10.0

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

n.d. n.d. n.d. 0.2 0.3 0.3 0.2 0.1 0.8 n.d. 0.6 0.6 n.d. n.d. n.d. n.d. 0.3 0.2 n.d. n.d. 0.1 n.d. n.d. n.d.

a All energies are in kcal mol−1, with error bars ±1 s.e.m. for calculated values. bBinding free energies from experimental instant inhibition assays ±1 s.d.52 cAverage binding free energies from preincubated assays ±1 s.e.m.47−51,62 n.d., not experimentally determined or reported.



annihilation of the additional methyl group in Ile. This was divided into three subperturbations, each of which consists of 50 evenly distributed FEP windows (λ-windows) of 30 ps: removal of partial charges, transformation of the regular van der Waals (Lennard−Jones) potential into a soft-core potential, and full annihilation of the atoms in the −CH3 group, followed by the creation of the corresponding hydrogen atom in valine.

RESULTS AND DISCUSSION

Conformational Analysis. The homology model of human COX-2 (hCOX-2) was built using the highest resolution crystal structure of a COX-2 enzyme, which corresponds to the murine species.25 (PDB code 3NT1, sequence identity = 87%). During the course of this project, crystal structures of hCOX-2 became available (PDB codes: 5F1A, 5IKQ, 5IKR, 5IKT, 5IKV, and 5KIR).44−46 Since the comparison of our homology model with

1240

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

Figure 4. (A) Plot of calculated vs experimental binding free energies of rapidly reversible (orange diamond), irreversible (green diamond), slowly reversible (blue diamonds), two-step slow tight-binding (purple diamonds), three-step slow tight-binding (red diamonds), other coxibs (red circles), and COX inhibitors with an unknown kinetic profile (gray circles). (B) Experimental54 and calculated relative binding free energies of the mutated hCOX-2 (V523I) compared to wt hCOX-2. Error bars show the calculated standard error of the mean.

diverse NSAIDs available, and we have previously demonstrated that including additional data points from other assays could be beneficial to cancel potential errors in the experimental results.22 Thus, where possible we use the average binding free energy from different sources (Tables 1 and 2), which include both preincubated48−51 and instant inhibition52 assays. The original data set contains 37 inhibitors with reported data on COX-2 inhibition. Each of them was docked in both the open and closed conformations and subjected to MD simulations to estimate the corresponding binding affinities using the LIE approach. Due to the size of the data set it would not be possible to discuss all the inhibitors in detail. Thus, we have chosen to focus on the data set of 13 of the inhibitors for which we know the kinetic properties (Table 1).15,53 However, computational results for the remaining 24 inhibitors are also given in Table 2. In this subset many, but not all, inhibitors converged to the closed conformation during MD simulations starting in the open conformation. The inhibitors that started and remained in the open conformation consistently had lower binding free energies in the open conformation than the closed one. Thus, it is likely that these inhibitors are time-dependent in COX-2. The reported binding free energies for this subset are extracted from the simulations that gave the lowest average binding free energies regardless of starting pose or protein conformation as stated in the column of preferred conformation (Table 2). The binding free energies from simulations in both conformations are reported for the other subset of 13 inhibitors, for which the kinetic class is known (Table 1). Figure 4A shows the correlation between the lowest binding free energies from MD simulations and experiment for all 37 inhibitors. As can be seen, the correlation is excellent with an average unsigned error of 0.44 kcal mol−1 for the smaller data set of 13 inhibitors, and 1.07 kcal mol−1 for the overall data set of 37 inhibitors, using the optimized γ = −4.9 (eq 8). We will proceed to discuss the binding modes and preferences of the smaller subset below. Free Energy Perturbation Calculations on the V523I Mutation in COX-2. Encouraged by the performance of our model we decided to perform in silico mutagenesis41−43 to test the hypothesis that the substitution of the gate-keeping residue

these structures gave an almost perfect superposition (RMSD of all heavy atoms between 0.65 and 0.73 Å), we decided to continue the ongoing molecular docking and MD simulations on our homology model. This model, as well as all crystal structures of human and murine COX-2, were in the closed conformation. The PMF calculations confirmed the relevance of these two conformations in COX-2, in analogy with experimental evidence18 and our previous calculations in COX-1.18,20,21 As with COX-1, the apo structure of COX-2 was found to be significantly more stable in the intrinsic closed conformation than in the open conformation (Figure 3A). However, there are important differences between the open conformations of COX-1 and COX-2. In COX-1, this conformation is characterized by the rotation of the Leu531 side-chain in the main binding pocket and a shift of helix D.18,21 The MD induced “open” conformation of COX-2 shows, in addition to these conformational changes described for COX-1, a backward displacement of Arg513, which thereby loses its coordination with Glu524 (Figure 3B). The latter residue retains a single salt bridge with Arg120 and creates a water channel that increases the solvation of the side pocket. This final conformation obtained from the PMF calculations is in agreement with the experimental mCOX-2 structure in complex with meloxicam, in the open conformation (PDB code 4M11, Figure S1).18,45 The PMF calculations of the rapidly reversible inhibitor ibuprofen in COX-2 reveal that the closed conformation remains the most stable one (Figure 3A), in agreement with our previous results in COX-1.20 A similar behavior was observed in the corresponding calculations in complex with the irreversible inhibitor acetylsalicylic acid. However, the open and closed conformations of COX-2 are equally stable in the presence of flurbiprofen (Figure 3A), a ligand that selectively stabilized the open conformation in COX-1.20,21 Thus, our calculations confirm that human COX-2 also has two stable conformations, where the open one is only induced by the binding of time-dependent inhibitors. Structure−Activity Relationships of NSAIDs. From the many experimental enzymatic assays reported for marketed NSAIDs, we chose to examine the data set of Warner and coworkers.47 This is the most extensive data set of structurally 1241

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

particular isoform, which agrees with the experimental lack of selectivity for this compound (Figure S2). Slow Irreversible Inhibitors. Aspirin is the only known truly irreversible inhibitor of COX in this data set. The acetylation by this compound of Ser530 has been well documented in both isoforms;17 thus the calculated binding affinities are not reflective of the irreversible step, but only of the initial binding step (Figure 6). Both our PMF (Figure 3A) and LIE calculations (Table 1) show that aspirin prefers binding to the closed conformation. Slow Reversible Inhibitors. The class of slow reversible inhibitors includes structurally diverse compounds. Some have the characteristic NSAID carboxylic acid moiety, while others are uncharged (Figure 2). This group of compounds is interesting because of the lower propensity of some inhibitors to cause gastrointestinal ulcers.56 In our simulations, we found that they display a variety of conformations. Some inhibitors, such as 6MNA and the structurally analogous naproxen (Figure 2) equally stabilize the open and closed conformations in regular MD simulations. Approximately equal numbers of simulations starting in the open conformation revert to the closed conformation, while simulations starting in the closed conformations convert into the open conformation. Similarly, meloxicam and piroxicam (Figure 2) are stable in both open and closed conformations, while others, such as the mefenamic acid and the COX-2 selective inhibitor nimesulide (Figure 2), do not stabilize the open conformation (Table 1). Thus, generalizations cannot be made for this class of inhibitors. However, we do find that for piroxicam, the open-closed affinities are very distinct and are well in agreement with the instant inhibition and preincubated assays (Table 1). Thus, the open-closed model remains applicable for the oxicams (Figure 7). Interestingly, nimesulide prefers the closed conformation, and it effectively destabilized the open conformation in all our simulations. Moreover, the calculated affinities in this closed conformation are well in agreement with observed affinities. Intriguingly, the corresponding binding free energies in the structurally similar NS-398 (Figure 2) correlated poorly with experimental binding free energies (Table 2), indicating that binding modes and affinities in this structural class might also be sensitive to small substitutions. Although both nimesulide and NS-398 are known to be time-dependent in hCOX-2 (but not in hCOX-1),15,57 it has not been established which of the kinetic classes NS-398 belongs to in hCOX-2. Experimentally, the V523I mutation nullifies the time-dependencies of both nimesulide and NS-398 and decreases their experimental affinities.58 The decrease in affinities is replicated in our V523I FEP calculations in complex with NS-398 (Figure 4B). However, at present we cannot explain the time-dependency of these two inhibitors solely using our open-closed model. Hence, the time-dependent property of these inhibitors requires some alternative explanation. To conclude, the reversibility seen in eq 3 could be explained by the open-closed model for some of the ligands in this class of inhibitors. In other ligands, the initial enzyme−inhibitor binding step may be time-dependent as proposed by others.59 In any case, this class of inhibitors needs further investigation to reveal the possible mechanisms and causes of time-dependency. Slow Tight-Binding Inhibitors (Two-Step Mechanism, eq 2). As opposed to the previously discussed kinetic classes, the slow tight-binding inhibitors show a clear distinction in affinities between the open and the closed conformations.

Ile523 for Val523 (Figure 1) would be the cause of selectivity between COX-1 and COX-2. This part of the study was done by FEP calculations where the wild-type COX-2 residue was mutated to the corresponding isoleucine (V523I mutant), in the presence and absence of four inhibitors belonging to different classes. The results, summarized in (Figure 4B), show that the inhibitor preference toward either amino acid is correctly captured by the FEP calculations. The resulting relative binding free energies showed that mefenamic acid (Figure 2) prefers the bulkier Ile523, while naproxen has a higher affinity with Val523, which agrees with the experimental mutational binding data.54 For both of these inhibitors the difference in binding free energy due to the mutation is relatively small (∼1 kcal mol−1). The calculations further show, in agreement with experiment, that binding of flurbiprofen is essentially unaffected by the mutation. Since this is known to be a preferential COX-1 inhibitor, the results clearly demonstrate that the gate-keeping residue in this case has no importance for selectivity. For the COX-2 selective inhibitor NS-398 the calculations predict a decrease in affinity due to the bulkier Ile523 side chain (Figure 4B). No inhibition of either the COX2 V523I mutant or native COX-1 was observed with NS-398 in ref 54, whereas other experiments on inhibition of COX-1 versus COX-2 indicated a 20-fold lower affinity of NS-398 for the latter enzyme,47,55 which is more in line with the above computational prediction. Rapidly Reversible Inhibitors. In this data set only one compound, ibuprofen (Figure 2), has been confirmed as a rapidly reversible inhibitor in COX-2. As expected, both the affinities and the kinetic profiles of this compound were consistent between the two COX isoforms (Figure S2). Ibuprofen prefers to bind in the native closed conformation (Figure 5) and, similar to our observations in COX-1, the

Figure 5. Average structure from MD simulations of the rapidly reversible inhibitor ibuprofen.

simulations starting in the open state rapidly converted to the closed conformation during the equilibration stage. This is in line with the PMF calculations, where the open conformation is destabilized compared to the corresponding closed conformation when the enzyme is in complex with ibuprofen (Figure 3A), again in agreement with the analogous calculations in COX-1.20 These data indicate ibuprofen should not prefer a 1242

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

Figure 6. Average structure from MD simulations of the proposed EI and EI* conformations of the irreversible inhibitor acetylsalicylic acid (aspirin).

Figure 7. Average structure from MD simulations of the proposed EI and EI* conformations of the slowly reversible inhibitor piroxicam with the crystal structure18 overlaid (gray).

both with regards to affinity and transition time, where COX-2 has a higher transition barrier and a lower affinity in the EI* conformation than COX-1 (eq 2).60 An explanation could be the increased space of the COX-2 binding pocket, where flurbiprofen can make translational movements toward the sidepocket (Figure 1) due to the substitution for the smaller Val523 in COX-2 compared to Ile523 in COX-1. As flurbiprofen enters the binding pocket the Leu531 side-chain rotates into the open conformation. This coincides with or precedes the influx of water molecules that line up between Tyr385 and Arg120 between Ser353 and the inhibitor. The waters form a stable water bridge that displaces flurbiprofen toward Leu531 and efficiently prevents the leucine side-chain from rotating back to the closed position (Figure 8A). The increased solvation of the binding pocket results in higher binding affinities of flurbiprofen in the open conformation compared to the closed one. This is similar, but not identical, to the mechanism in COX-1 where a water bridge forms between Ser530 and Arg120 on the opposite side of the inhibitor.20,21 However, due to the reduced tightness of the binding pocket in COX-2, and the resulting mobility of flurbiprofen, the water bridge in COX2 is less stable than in COX-1.21 With the reduced stability of the water bridge Leu531 can more readily rotate back into the closed conformation. Thus, the EI → EI* transition, or at least the rotational state of Leu531, is more reversible in COX-2 than in COX-1, affecting the observed affinity in COX-2.

The open conformation was the lowest energy conformation for all two-step slow tight-binding inhibitors (eq 4). These were in agreement with the binding affinities of preincubated assays, while the binding free energies of the closed conformation simulations reproduced the binding affinities of instantaneous inhibition assays (Table 1). This is in line with previous results in COX-1.20 From these calculations we can clearly see that, similar to hCOX-1, the affinities of these inhibitors favor the open conformation of hCOX-2. However, for all inhibitors in Table 1, whenever the open conformation yields the lower binding free the affinity increase in the open conformation is smaller in COX-2 than in COX-1 for all probed inhibitors (Table 1, Figure S2). For instance, our LIE calculations provided an average binding free energy of flurbiprofen (Figure 2) in COX-2 of 8.8 kcal·mol−1 in the open and 7.8 kcal·mol−1 in the closed conformation (8.2 and 7.8 experimentally). The corresponding calculation in COX-1 yielded 10.6 kcal·mol−1 in the open and 9.1 kcal·mol−1 in the closed conformations.20 Although flurbiprofen is the inhibitor that best stabilizes the open conformation in PMF calculations (Figure 3A), it is only able to stabilize the open conformation to approximately the same degree as the closed conformation. This is in contrast to in COX-1, where the open conformation was stabilized more than the closed conformation.20,21 The relative difference is ∼2.5 kcal·mol−1 between the two open conformations (Figure S2). Experimentally, a difference has been repeatedly observed 1243

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

Figure 8. Average structures from MD simulations of the proposed EI and EI* conformations of two-step slow tight-binding inhibitors: (A) flurbiprofen, (B) diclofenac.

The preferential binding toward either COX-1 or COX-2 in this class of ligands depends on the ability of the inhibitor to use the increased space created in the binding pocket in the open conformation. Flurbiprofen, which relies on the previously described stabilized water bridge to retain the open conformation, has less affinity for COX-2 than COX-1. The substitution of Ile523 for Val523 in COX-2 results in increased space in COX-2 (Figure 1), which makes the inhibitor more mobile, leading to a less stable water bridge, a less stabilized open conformation, and consequently lower binding affinities. In this complex, we also probed the effect of mutating the Val523 to see how the point mutation would affect binding (Figure 4B). Surprisingly, the effect is negligible, which suggests that the main effect of selectivity in flurbiprofen is due to the open and closed conformation rather than the V523I substitution. From our MD simulations exclusively in the open conformation, we obtain slightly higher affinities (0.6 kcal mol−1) than indicated by the experimental data. As noted before, occasional rotations of this isoleucine after the other conformational changes have been established do not lead to the total recovery of the intrinsic closed conformation. Thus, under physiological conditions mixed rotational conformations of Leu531 are probably present. Diclofenac, on the other hand, is bulkier and more flexible than flurbiprofen. It benefits more from the increased solvation and space in the larger side pocket that is created through the substitution of Ile523 in COX-1 to Val523 in COX-2. Slow Tight-Binding Inhibitor (Three-Step Mechanism, eq 4). In this data set only one inhibitor, celecoxib (Figure 2), has

Although most two-step slow tight-binding inhibitors bind preferentially to COX-1, there is one exception: diclofenac (Figure 8B). Experimentally, diclofenac has ∼0.8 kcal mol−1 higher affinity toward hCOX-2 than hCOX-1 in preincubated assays (Figure S2). This difference is due to the increased space in the gate-area between the main and side-pocket, regulated by V523 in COX-2 and I523 in COX-1. The smaller residue allows the bulky diclofenac to position itself optimally toward the Arg120 side-chain (Figure 8B). This result agrees with experimental studies of the V523I mutation, which resulted in an almost 1 kcal mol−1 decrease in affinity for diclofenac.54 Apart from the differences between the two isoforms, diclofenac also shows a marked difference in binding free energies between the open and closed conformation (4.2 kcal mol−1), in agreement with experimental results of preincubated and instant inhibition assays, respectively (Table 1). The energetic difference between the two conformations is larger in COX-2 than in COX-1, and originates predominantly from more favorable electrostatics in the open conformation, while the nonpolar interactions differ by less than 1 kcal/mol. Besides enhanced interactions with Arg120, the influx of water molecules increases solvation of the otherwise hydrophobic binding pocket. Notably, a stable water bridge is formed in the open conformation, similar to the previously discussed case of flurbiprofen, but stable water molecules are also present between the diclofenac carboxylate and Tyr355. The same behavior is also seen for the two-step slow tight-binding inhibitor indomethacin, for which the binding free energies increase by 3.8 kcal mol−1 in the open conformation. 1244

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

Figure 9. Average structures from MD simulations of the proposed EI, EI*, and EX conformations of the three-step slow tight-binding inhibitor celecoxib.

the demonstrated kinetic profile of a three-step slow tightbinding inhibitor. Celecoxib has previously been shown to bind in the side pocket of COX in the closed conformation (PDB codes 3KK6 and 3LN1)14,61 Here, we simulated it bound in the side pocket in both the open and the closed conformation. Two binding modes of celecoxib in the main binding pocket starting in both open and closed conformations, but without reaching into the side pocket were also probed. In the main pocket, all simulations of coxibs converged into the closed conformation. Furthermore, the conformation in which the fluorines were oriented toward Arg120 had the highest affinities. In the side pocket, the inhibitors were docked as in reported crystal structures.14,61 In the closed conformation, the inhibitor remained mostly in the initial pose, but it shifted slightly in the open conformation and displayed a novel conformation (Figure 9) during the equilibration stage of unbiased MD simulations. This conformation was very stable, and as the enzyme stayed in the open conformation, the inhibitor remained there for the duration of the 10 ns of MD simulations, with the resulting binding free energies well in agreement with experimental data (Table 3). The main difference between the two complexes is found in the enzyme−ligand electrostatic interactions. In the open conformation the celecoxib sulfone oxygen can form hydrogen bonds with both the Arg513 side-chain and the Tyr355 hydroxyl group, while the nitrogen interacts with Glu524. In the closed conformation only one hydrogen bond is formed

Table 3. Average Energies from MD Simulations in Different Binding Modesa ligand celecoxib

conformation water main pocket (closed) side pocket (closed) side pocket (open)

⟨Uvdw l‑s ⟩

⟨Uell‑s⟩

ΔGcalc bind

−27.7 ± 0.0 −44.8 ± 0.9

−43.6 ± 0.1 −33.4 ± 1.1

−3.5 ± 0.5

−48.6 ± 1.2

−33.6 ± 1.0

−4.3 ± 0.5

−46.6 ± 0.2

−44.1 ± 0.3

−9.2 ± 0.1

a All energies are in kcal mol−1, with error bars ±1 s.e.m. for calculated values.

(Figure 9), thus yielding poorer electrostatic energies (Table 3). Given that Arg513 is not present in COX-1 (where it is substituted for His513), selectivity of the coxibs toward COX-2 could be explained by a combination of the wider entrance caused by V523I substitution, and the optimal interactions in the open conformation with the displacement of Arg513, which cannot occur in COX-1. The isoleucine in COX-1 creates a smaller gate to the side pocket, which would reduce the flexibility and mobility of the inhibitor. In combination with the reduced flexibility of His513, this hampers the possibility of obtaining optimal interactions with His513. Thus, we find that it is a combination of the time-dependent open conformation, and the existence of Arg513 and Val523 in the side pocket that 1245

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Biochemistry



creates the optimal interactions for coxibs, thus resulting in selectivity toward COX-2. Kinetic experiments for celecoxib indicate that the timedependent IC50 corresponds to a ∼5 kcal mol−1 higher affinity than the initial inhibitory interaction.53 This is in agreement with our calculated values, where the binding in the main −1 pocket yields a difference in ΔGcalc bind of ∼5.7 kcal mol . The difference in affinity with the initial binding in the closed conformation side pocket yields ∼4.9 kcal mol−1. Similar poses were then tried for the two other coxibs known to bind in the side pocket (Figure 1), rofecoxib and sc-58125 (Figure 2, Table 2). The results were similar to those with celecoxib, with around 4 kcal mol−1 increase in affinities in the open conformation due to the displacement of Arg513. Thus, we can propose a model in which the binding in the main binding pocket corresponds to the EI conformation, the initial binding in the side pocket corresponds to the EI* conformation, and the binding in the open conformation side pocket corresponds to the EX conformation of eq 4 (Figure 9).

CONCLUSION Using various computational and free energy methods we find that the conformational change regulating the two-step, slow tight-binding inhibition in COX-1 can also be induced in COX2 for the same class of inhibitors. The same conformational change is seen in COX-2 selective coxibs following the threestep reaction mechanism and for some, but not all, inhibitors in the class of slow reversible inhibitors. Thus, MD simulations and free energy calculations support a three-step mechanism for COX-2 selective coxibs where the initial binding occurs in the main pocket. The subsequent reversible insertion into the side pocket is then followed by the slowly time-dependent binding in the side pocket where Arg513 has been displaced. ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.biochem.7b01019. Figures comparing open and closed crystal conformations and experimental binding free energies of inhibitors in COX-1 and COX-2 (PDF)



REFERENCES

(1) Picot, D., Loll, P. J., and Garavito, R. M. (1994) The X-ray crystal structure of the membrane protein prostaglandin H2 synthase-1. Nature 367, 243−249. (2) Kurumbail, R. G., Stevens, A. M., Gierse, J. K., and McDonald, J. J. (1996) Structural basis for selective inhibition of cyclooxygenase-2 by anti-inflammatory agents. Nature 384, 644−648. (3) Copeland, R. A., Williams, J. M., Giannaras, J., Nurnberg, S., Covington, M., Pinto, D., Pick, S., and Trzaskos, J. M. (1994) Mechanism of selective inhibition of the inducible isoform of prostaglandin G/H synthase. Proc. Natl. Acad. Sci. U. S. A. 91, 11202−11206. (4) Steinmeyer, J. (2000) Pharmacological basis for the therapy of pain and inflammation with nonsteroidal anti-inflammatory drugs. Arthritis Res. 2, 379−385. (5) Winkelmayer, W. C., Waikar, S. S., Mogun, H., and Solomon, D. H. (2008) Nonselective and cyclooxygenase-2-selective NSAIDs and acute kidney injury. Am. J. Med. 121, 1092−1098. (6) Schneider, V., Lévesque, L. E., Zhang, B., Hutchinson, T., and Brophy, J. M. (2006) Association of selective and conventional nonsteroidal antiinflammatory drugs with acute renal failure: A population-based, nested case-control analysis. Am. J. Epidemiol. 164, 881−889. (7) Bresalier, R. S., Sandler, R. S., Quan, H., Bolognese, J. A., Oxenius, B., Horgan, K., Lines, C., Riddell, R., Morton, D., Lanas, A., Konstam, M. A., and Baron, J. A. (2006) Cardiovascular events associated with rofecoxib in a colorectal adenoma chemoprevention trial. N. Engl. J. Med. 355, 221−221. (8) Marnett, L. J., and DuBois, R. N. (2002) COX-2: A target for colon cancer prevention. Annu. Rev. Pharmacol. Toxicol. 42, 55−80. (9) Brown, J. R., and DuBois, R. N. (2005) COX-2: a molecular target for colorectal cancer prevention. J. Clin. Oncol. 23, 2840−2855. (10) Soslow, R. A., Dannenberg, A. J., Rush, D., Woerner, B. M., Khan, K. N., Masferrer, J., and Koki, A. T. (2000) COX-2 is expressed in human pulmonary, colonic, and mammary tumors. Cancer 89, 2637−2645. (11) Tucker, O. N., Dannenberg, A. J., Yang, E. K., Zhang, F., Teng, L., Daly, J. M., Soslow, R. A., Masferrer, J. L., Woerner, B. M., Koki, A. T., and Fahey, T. J. (1999) Cyclooxygenase-2 expression is upregulated in human pancreatic cancer. Cancer Res. 59, 987−990. (12) Sano, H., Kawahito, Y., Wilder, R. L., Hashiramoto, A., Mukai, S., Asai, K., Kimura, S., Kato, H., Kondo, M., and Hla, T. (1995) Expression of cyclooxygenase-1 and −2 in human colorectal cancer. Cancer Res. 55, 3785−3789. (13) Hawkey, C. J. (1999) COX-2 inhibitors. Lancet 353, 307−314. (14) Rimon, G., Sidhu, R. S., Lauver, D. A., Lee, J. Y., Sharma, N. P., Yuan, C., Frieler, R. A., Trievel, R. C., Lucchesi, B. R., and Smith, W. L. (2010) Coxibs interfere with the action of aspirin by binding tightly to one monomer of cyclooxygenase-1. Proc. Natl. Acad. Sci. U. S. A. 107, 28−33. (15) Walker, M. C., Kurumbail, R. G., Kiefer, J. R., Moreland, K. T., Koboldt, C. M., Isakson, P. C., Seibert, K., and Gierse, J. K. (2001) A three-step kinetic mechanism for selective inhibition of cyclooxygenase-2 by diarylheterocyclic inhibitors. Biochem. J. 357, 709−718. (16) So, O.-Y., Scarafia, L. E., Mak, A. Y., Callan, O. H., and Swinney, D. C. (1998) The Dynamics of Prostaglandin H Synthases. J. Biol. Chem. 273, 5801−5807. (17) Lei, J., Zhou, Y., Xie, D., and Zhang, Y. (2015) Mechanistic insights into a classic wonder drug–aspirin. J. Am. Chem. Soc. 137, 70− 73. (18) Xu, S., Hermanson, D. J., Banerjee, S., Ghebreselasie, K., Clayton, G. M., Garavito, R. M., and Marnett, L. J. (2014) Oxicams bind in a novel mode to the cyclooxygenase active site via a two-watermediated h-bonding network. J. Biol. Chem. 289, 6799−6808. (19) Vecchio, A. J., Simmons, D. M., and Malkowski, M. G. (2010) Structural basis of fatty acid substrate binding to cyclooxygenase-2. J. Biol. Chem. 285, 22152−22163.





Article

AUTHOR INFORMATION

Corresponding Author

*Phone: +46 18 471 4109. Fax: +46 18 53 69 71. E-mail: [email protected]. ORCID

Hugo Gutiérrez-de-Terán: 0000-0003-0459-3491 Johan Åqvist: 0000-0003-2091-0610 Funding

Support from the Swedish Research Council, the eSSENCE escience initiative, and the Swedish National Infrastructure for Computing (SNIC) is gratefully acknowledged. Notes

The authors declare no competing financial interest.



ABBREVIATIONS NSAID, nonsteroidal anti-inflammatory drug; MD, molecular dynamics; FEP, free energy perturbation; LIE, linear interaction energy; COX, cyclooxygenase; PMF, potential of mean force 1246

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry

(41) Boukharta, L., Gutiérrez-de-Terán, H., and Åqvist, J. (2014) Computational prediction of alanine scanning and ligand binding energetics in G-protein coupled receptors. PLoS Comput. Biol. 10, e1003585. (42) Keränen, H., Gutiérrez-de-Terán, H., and Åqvist, J. (2014) Structural and energetic effects of A2A adenosine receptor mutations on agonist and antagonist binding. PLoS One 9, e108492. (43) Keränen, H., Åqvist, J., and Gutiérrez-de-Terán, H. (2015) Free energy calculations of A2A adenosine receptor mutation effects on agonist binding. Chem. Commun. 51, 3522−3525. (44) Lucido, M. J., Orlando, B. J., Vecchio, A. J., and Malkowski, M. G. (2016) Crystal structure of aspirin-acetylated human cyclooxygenase-2: Insight into the formation of products with reversed stereochemistry. Biochemistry 55, 1226−1238. (45) Orlando, B. J., and Malkowski, M. G. (2016) Crystal structure of rofecoxib bound to human cyclooxygenase-2. Acta Crystallogr., Sect. F: Struct. Biol. Commun. 72, 772−776. (46) Orlando, B. J., and Malkowski, M. G. (2016) Substrate-selective inhibition of cyclooxygeanse-2 by fenamic acid derivatives is dependent on peroxide tone. J. Biol. Chem. 291, 15069−15081. (47) Warner, T. D., Giuliano, F., Vojnovic, I., Bukasa, A., Mitchell, J. A., and Vane, J. R. (1999) Nonsteroid drug selectivities for cyclooxygenase-1 rather than cyclo-oxygenase-2 are associated with human gastrointestinal toxicity: A full in vitro analysis. Proc. Natl. Acad. Sci. U. S. A. 96, 7563−7568. (48) O’Neill, G. P., Mancini, J. A., Kargman, S., Yergey, J., Kwan, M. Y., Falgueyret, J. P., Abramovitz, M., Kennedy, B. P., Ouellet, M., and Cromlish, W. (1993) Overexpression of human prostaglandin G/H synthase-1 and -2 by recombinant vaccinia virus: inhibition by nonsteroidal anti-inflammatory drugs and biosynthesis of 15hydroxyeicosatetraenoic acid. Mol. Pharmacol. 45, 245−254. (49) Glaser, K., Sung, M.-L., O’Neill, K., Belfast, M., Hartman, D., Carlson, R., Kreft, A., Kubrak, D., Hsiao, C.-L., and Weichman, B. (1995) Etodolac selectively inhibits human prostaglandin G/H synthase 2 (PGHS-2) versus human PGHS-1. Eur. J. Pharmacol. 281, 107−111. (50) Brideau, C., Kargman, S., Liu, S., Dallob, A. L., Ehrich, E. W., Rodger, I. W., and Chan, C. C. (1996) A human whole blood assay for clinical evaluation of biochemical efficacy of cyclooxygenase inhibitors. Inflammation Res. 45, 68−74. (51) Cryer, B., and Feldman, M. (1998) Cyclooxygenase-1 and Cyclooxygenase-2 Selectivity of Widely Used Nonsteroidal AntiInflammatory Drugs. Am. J. Med. 104, 413−421. (52) Laneuville, O., Breuer, D. K., Dewitt, D. L., Hla, T., Funk, C. D., and Smith, W. L. (1994) Differential inhibition of human prostaglandin endoperoxide H synthases-1 and-2 by nonsteroidal anti-inflammatory drugs. J. Pharmacol. Exp. Ther. 271, 927−934. (53) Gierse, J. K., Koboldt, C. M., Walker, M. C., Seibert, K., and Isakson, P. C. (1999) Kinetic basis for selective inhibition of cyclooxygenases. Biochem. J. 339, 607−614. (54) Gierse, J. K., McDonald, J. J., Hauser, S. D., Rangwala, S. H., Koboldt, C. M., and Seibert, K. (1996) A single amino acid difference between cyclooxygenase-1 (COX-1) and −2 (COX-2) reverses the selectivity of COX-2 specific inhibitors. J. Biol. Chem. 271, 15810− 15814. (55) Kato, M., Nishida, S., Kitasato, H., Sakata, N., and Kawai, S. (2001) Cyclooxygenase-1 and cyclooxygenase-2 selectivity of nonsteroidal anti-inflammatory drugs: investigation using human peripheral monocytes. J. Pharm. Pharmacol. 53, 1679−1685. (56) Shah, A. A., Thjodleifsson, B., Murray, F. E., Kay, E., Barry, M., Sigthorsson, G., Gudjonsson, H., Oddsson, E., Price, A. B., Fitzgerald, D. J., and Bjarnason, I. (2001) Selective inhibition of COX-2 in humans is associated with less gastrointestinal injury: a comparison of nimesulide and naproxen. Gut 48, 339−346. (57) Ouellet, M., and Percival, M. D. (1995) Effect of inhibitor timedependency on selectivity towards cyclooxygenase isoforms. Biochem. J. 306 (Pt 1), 247−251. (58) Guo, Q., Wang, L.-H., Ruan, K.-H., and Kulmacz, R. J. (1996) Role of Val509 in time-dependent inhibition of human prostaglandin

(20) Shamsudin Khan, Y., Kazemi, M., Gutiérrez-de-Terán, H., and Åqvist, J. (2015) Origin of the enigmatic stepwise tight-binding inhibition of cyclooxygenase-1. Biochemistry 54, 7283−7291. (21) Shamsudin Khan, Y., Gutiérrez-de-Terán, H., and Åqvist, J. (2017) Probing the time dependency of cyclooxygenase-1 inhibitors by computer simulations. Biochemistry 56, 1911−1920. (22) Shamsudin Khan, Y., Gutiérrez-de-Terán, H., Boukharta, L., and Åqvist, J. (2014) Toward an optimal docking and free energy calculation scheme in ligand design with application to COX-1 inhibitors. J. Chem. Inf. Model. 54, 1488−1499. (23) Eswar, N., Webb, B., Marti-Renom, M. A., Madhusudhan, M. S., Eramian, D., Shen, M.-Y., Pieper, U., and Sali, A. (2006) Comparative protein structure modeling using Modeller. Curr. Protoc Bioinformatics, 5.6.1. (24) Gupta, K., Selinsky, B. S., Kaub, C. J., Katz, A. K., and Loll, P. J. (2004) The 2.0 Å resolution crystal structure of prostaglandin H2 synthase-1: Structural insights into an unusual peroxidase. J. Mol. Biol. 335, 503−518. (25) Duggan, K. C., Walters, M. J., Musee, J., Harp, J. M., Kiefer, J. R., Oates, J. A., and Marnett, L. J. (2010) Molecular basis for cyclooxygenase inhibition by the non-steroidal anti-inflammatory drug naproxen. J. Biol. Chem. 285, 34950−34959. (26) Schrödinger Release 2011-2: MacroModel; Schrödinger, LLC: New York, 2011. (27) Schrödinger Release 2011-2: Maestro; Schrödinger, LLC: New York, 2011. (28) Friesner, R. A., Murphy, R. B., Repasky, M. P., Frye, L. L., Greenwood, J. R., Halgren, T. A., Sanschagrin, P. C., and Mainz, D. T. (2006) Extra precision glide: Docking and scoring incorporating a model of hydrophobic enclosure for protein−ligand complexes. J. Med. Chem. 49, 6177−6196. (29) Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983) Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926−935. (30) Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J. (1996) Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225−11236. (31) Marelius, J., Kolmodin, K., Feierberg, I., and Åqvist, J. (1998) Q: a molecular dynamics program for free energy calculations and empirical valence bond simulations in biomolecular systems. J. Mol. Graphics Modell. 16, 213−225. (32) King, G., and Warshel, A. (1989) J. Chem. Phys. 91, 3647−3661. (33) Ryckaert, J.-P., Ciccotti, G., and Berendsen, H. J. C. (1977) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327−341. (34) Lee, F. S., and Warshel, A. (1992) A local reaction field method for fast evaluation of long-range electrostatic interactions in molecular simulations. J. Chem. Phys. 97, 3100−3107. (35) Wallin, G., Härd, T., and Åqvist, J. (2012) Folding-reaction coupling in a self-cleaving protein. J. Chem. Theory Comput. 8, 3871− 3879. (36) Åqvist, J., Medina, C., and Samuelsson, J.-E. (1994) A new method for predicting binding affinity in computer-aided drug design. Protein Eng., Des. Sel. 7, 385−391. (37) Hansson, T., Marelius, J., and Åqvist, J. (1998) Ligand binding affinity prediction by linear interaction energy methods. J. Comput.Aided Mol. Des. 12, 27−35. (38) Almlöf, M., Carlsson, J., and Åqvist, J. (2007) Improving the accuracy of the linear interaction energy method for solvation free energies. J. Chem. Theory Comput. 3, 2162−2175. (39) Marelius, J., Hansson, T., and qvist, J. (1998) Calculation of ligand binding free energies from molecular dynamics simulations. Int. J. Quantum Chem. 69, 77−88. (40) Cheng, Y.-C., and Prusoff, W. H. (1973) Relationship between the inhibition constant (KI) and the concentration of inhibitor which causes 50% inhibition (I50) of an enzymatic reaction. Biochem. Pharmacol. 22, 3099−3108. 1247

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248

Article

Biochemistry H synthase-2 cyclooxygenase activity by isoform-selective agents. J. Biol. Chem. 271, 19134−19139. (59) Blobaum, A. L., Xu, S., Rowlinson, S. W., Duggan, K. C., Banerjee, S., Kudalkar, S. N., Birmingham, W. R., Ghebreselasie, K., and Marnett, L. J. (2015) Action at a distance: Mutations of peripheral residues transform rapid reversible inhibitors to slow, tight binders of cyclooxygenase-2. J. Biol. Chem. 290, 12793−12803. (60) Callan, O. H., So, O. Y., and Swinney, D. C. (1996) The kinetic factors that determine the affinity and selectivity for slow binding inhibition of human prostaglandin H synthase 1 and 2 by indomethacin and flurbiprofen. J. Biol. Chem. 271, 3548−3554. (61) Wang, J. L., Limburg, D., Graneto, M. J., Springer, J., Hamper, J. R. B., Liao, S., Pawlitz, J. L., Kurumbail, R. G., Maziasz, T., Talley, J. J., Kiefer, J. R., and Carter, J. (2010) The novel benzopyran class of selective cyclooxygenase-2 inhibitors. Part 2: The second clinical candidate having a shorter and favorable human half-life. Bioorg. Med. Chem. Lett. 20, 7159−7163. (62) Houtzager, V., Ouellet, M., Falgueyret, J. P., Passmore, L. A., Bayly, C., and Percival, M. D. (1996) Inhibitor-induced changes in the intrinsic fluorescence of human cyclooxygenase-2. Biochemistry 35, 10974−10984.

1248

DOI: 10.1021/acs.biochem.7b01019 Biochemistry 2018, 57, 1236−1248