Atomistic Insights into the Inhibition of Cysteine Proteases: First QM

Apr 8, 2008 - Epoxides and aziridines are important building blocks for inhibitors of cysteine proteases which are promising drug targets for many dis...
0 downloads 0 Views 625KB Size
Download PDF
5458

J. Phys. Chem. B 2008, 112, 5458-5469

Atomistic Insights into the Inhibition of Cysteine Proteases: First QM/MM Calculations Clarifying the Regiospecificity and the Inhibition Potency of Epoxide- and Aziridine-Based Inhibitors Milena Mladenovic,† Konstantin Junold,† Reinhold F. Fink,† Walter Thiel,§ Tanja Schirmeister,‡ and Bernd Engels*,† Institut fu¨r Organische Chemie, and Institut fu¨r Pharmazie und Lebensmittelchemie, UniVersita¨t Wu¨rzburg, Am Hubland, D-97074 Wu¨rzburg, Germany, and Max-Planck-Institut fu¨r Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mu¨lheim an der Ruhr, Germany ReceiVed: NoVember 29, 2007; In Final Form: January 23, 2008

Epoxides and aziridines are important building blocks for inhibitors of cysteine proteases which are promising drug targets for many diseases. In spite of the large amount of experimental data concerning inhibition potency, structure-activity relationships, and structural arrangements of enzyme-inhibitor complexes, little is known about the basic principles which connect the substitution pattern with the resulting activities. To shed some light on this issue which is essential for the rational design of improved compounds, we have studied the inhibition processes theoretically for various inhibitors using quantum mechanical/molecular mechanical hybrid approaches and classical molecular dynamics simulations. The careful analysis of the computational results allows insight into the interactions which govern the regio- and stereospecificity of the interactions. Known structure-activity relationships are rationalized in terms of the same interactions that determine the measured pH dependencies. Inconsistencies in existing X-ray structures are resolved through comparison with the computed structures, which leads to a reassessment of the factors that control the inhibition potency. Similarities and differences in the mode of action of epoxide- and aziridine-based inhibitors are elucidated. Finally the small reaction barriers computed for the irreversible step in E64 analogues call into question the commonly accepted two-step model of inhibition since the second, irreversible step is predicted to be so fast that suitably oriented enzyme-inhibitor complexes will react rather than dissociate and equilibrate.

1. Introduction Epoxides and aziridines are important in polymer,1 supramolecular,2 and organic chemistry.3,4 In medicinal chemistry epoxysuccinic acid derivatives and analogous aziridines serve as building blocks for inhibitors of cysteine proteases,5 which represent promising drug targets for many diseases.6-12 Some epoxysuccinic acid derivatives have even reached clinical development as anticancer drugs.13 Aziridine-2,3-dicarboxylic acid derivatives have proven to be selective, nontoxic cysteine protease inhibitors suitable as antiparasitic14 and affinity-labeling agents.15 Containing a structural motif found in a natural compound (E-64, Figure 1), already in the late 1970s,16 these inhibitor classes have gained much recent interest as documented by more than 120 publications and patents during the last 3 years. Epoxides and aziridines are irreversible inhibitors of cysteine proteases. Their mode of action is commonly represented by a two-step reaction: Ki

ki

E + I {\} EI 98 E-I In the first step, characterized by the dissociation constant Ki, a reversible noncovalent enzyme inhibitor complex EI is formed. * Corresponding author. E-mail: [email protected]. Phone: ++49(0)9318885394. † Institut fu ¨ r Organische Chemie, Universita¨t Wu¨rzburg. ‡ Institut fu ¨ r Pharmazie und Lebensmittelchemie, Universita¨t Wu¨rzburg. § Max-Planck-Institut fu ¨ r Kohlenforschung.

The second step, characterized by the first-order rate constant of inhibition ki, is initiated by the nucleophilic attack of the negatively charged Cys thiolate at one of the electrophilic ring carbons (Figure 1). The ring-opening reaction leads to irreversible alkylation of the cysteine protease. The second-order rate constant k2nd ) ki/Ki is a measure for the overall inhibition potency of such inhibitors. According to the X-ray structures of the products E-I with epoxysuccinyl-based inhibitors17-22 and to model reactions with aziridine-2,3-dicarboxylates,23 the ring opening proceeds as an SN2 reaction. The overall inhibition process is determined by a complicated interplay between substituents and environmental effects. The experimental findings can be summarized as follows: (i) The known chemical difference between epoxides and aziridines in reactivity against nucleophilic ring opening is reflected in their inhibition potency: for pH values around 7 the aziridines are less effective, but for pH values around 4, N-unsubstituted and N-alkylated aziridines are as potent as epoxides. The lower inhibition potency at pH 6-7 is mainly due to lower ki values. (ii) For both epoxide- and aziridinebased inhibitors a free acid attached to the three-membered ring enhances inhibition potency compared to esters.24-26 However, this is only the case for compounds which bind solely into the nonprimed binding sites of the proteases (E-64 analogues). (iii) The regioselectivity of the alkylation step has been elucidated by X-ray investigations with epoxide-based inhibitors: in compounds bearing an acid at C2 only the R-carbon (C2) is attacked (Figure 1),17-20 whereas in the corresponding esters

10.1021/jp711287c CCC: $40.75 © 2008 American Chemical Society Published on Web 04/08/2008

QM/MM Calculations of Epoxide-Based Inhibitors

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5459

Figure 1. Left: inactivation of cysteine proteases by epoxysuccinyl peptides and corresponding aziridines with free acid at C2. The atom numbering corresponds to the 1ITO pdb structure (ref 17). The QM part is indicated in bold. Right: ring opening of R,β-epoxycarbonyl compounds by methyl thiolate in solution (ref 29).

and amides both attack at C2 or C3 is observed, depending on the overall binding mode.27 (iv) An experimental differentiation between the reversible (Ki) and irreversible (ki) step of the inhibition has been achieved by Meara and Rich with E-64 analogues.24,28 The following order of inhibition activity concerning the C2 substituent is found:28 CO2H > CONH2 > COCH3 > CO2Et ≈ CO2Me > CH2OH ≈ H. The replacement of the acid (CO2H) by hydroxamic acid (CONHOH), amide (CONH2), or ketone (COCH3) mainly influences Ki, whereas ki remains nearly unaffected. (v) These results for the ring opening of epoxide-based inhibitors by the Cys residue of the proteases are not in line with those for the ring opening of R,βepoxycarbonyl compounds by methyl thiolate in solution (Figure 1, right-hand side).29 For these latter reactions Bihovsky observed the following trend for the rate constant k: COCH3 > CO2CH3 > CONH2 > CO2-. These investigations also revealed a different regioselectivity. For esters the rate constants for R- and β-attack are quite similar (kR ) 0.25 min-1 vs kβ ) 0.19 min-1), whereas for acids they differ by a factor of 2 (kR ) 0.0036 min-1 vs kβ ) 0.0064 min-1) with the β-attack preferred. In spite of this large body of experimental data, much is still unknown about the basic principles which govern the inhibition potencies and the structure-activity relationships (SAR). Since the enzyme environment is expected to possess a strong influence on the inhibition process it must be taken into account to obtain a realistic theoretical description. However, although computations about the mechanism of the hydrolysis of a peptide were performed,30 corresponding investigations about possible inhibition mechanisms are missing. At hand are only quantum mechanical calculations23,31,32 which have provided a detailed picture for the ring-opening reactions with methyl thiolate in a polar solvent. In the present paper these enzymatic reactions are studied by means of hybrid methods33 (quantum mechanics/ molecular modeling, QM/MM) that combine quantum mechanical methods (QM) for the description of bond-forming and bond-breaking processes and molecular mechanics approaches (MM) for the consideration of the influence of the protein environment and the surrounding solvent. In addition, dynamical effects are investigated by classical molecular dynamics (MD) simulations. By comparing various model systems the calculations allow a detailed analysis of kinetics and thermodynamics of the mode of action of epoxide- and aziridine-based inhibitors and elucidate the main interactions governing these processes. We address regioselectivity and stereoselectivity as well as SAR and the pH dependency of the inhibition. Our calculations

indicate inconsistencies in the previous X-ray structures and question the generally accepted two-step model of inhibition. Furthermore, they rationalize the similarities and differences in the mode of action of epoxide- and aziridine-based inhibitors. The article is organized as followed: In section 2, we briefly describe the chosen computational methods and model systems. In section 3, we present and interpret the computational results for the epoxide-based and aziridine-based inhibitors. Section 4 offers a summary of the major findings. 2. Computational Methods and Model Systems The irreversible step in the inhibition reaction of cysteine proteases takes place in the active site of the enzyme and involves the creation of the covalent bond between the sulfur atom of the Cys29 residue and the carbon atom of the inhibitor ring, as well as breaking of a C-O or C-N bond of the threemembered ring during the ring opening. Additionally, a protonation of the emerging alcoholate or amide moiety can occur. The enzymatic environment plays also an important role in the inhibition, especially the nearby oxyanion hole (side chain of Gln23 and backbone NH group of Cys29, Figure 1) which stabilizes the enzyme-inhibitor complex by noncovalent interactions. Finally, the solvent has to be taken into account since the active site lies on the surface of the enzyme and is solvent accessible. Since such large systems cannot be treated quantum mechanically, the method of choice is the combined QM/MM approach. QM/MM methods34 divide the total system (enzyme, solvent, and inhibitor) into the active center and the rest. The active center comprises all parts which are directly involved in the bond-breaking/formation processes and thus have to be described by QM approaches, whereas the influence of the protein environment and the solvent can be captured at the MM level. The QM and MM regions interact by electrostatic and van der Waals terms. In our work the electrostatic QM/MM interactions are represented by an electronic embedding scheme35 which incorporates the MM charges into the one-electron QM Hamiltonian and thereby allows the QM system to adapt to the field exerted by the environment. Dangling bonds at the QM/MM boundary are capped with hydrogen link atoms36 in the framework of the charge shift method. The overall system used as a model to study the covalent bond formation in the cysteine protease inhibition reaction by epoxide- or aziridine-based inhibitors comprises the cysteine protease, the inhibitor, and additional water molecules. Starting

5460 J. Phys. Chem. B, Vol. 112, No. 17, 2008 structures for the computations were obtained by the following preparation procedure. In the X-ray structure of the cathepsin B-E64c complex (1ITO)17 the E64c was replaced by the respective inhibitor model (see below). This system was solvated in a water sphere with the radius of 50 Å. The charge of the resulting structure was compensated by seven sodium atoms which were placed at the outer surface of the solvent shell. The whole system was then submitted to a series of consecutive constrained optimizations, MD runs (heating from 100 to 300 K and equilibration at 300 K; time step 1 fs), and further solvation until saturation with regard to the number of water molecules was reached. In this procedure the enzyme was allowed to “move” only during the optimization and was held fixed during the MD runs. The outer 5 Å layer of water and the inhibitor were kept fixed all the time. The preparation of the system was performed with the CHARMM program package.37 Water molecules were of TIP3 kind. In the subsequent QM/MM calculations only the inner water sphere of r ) 25 Å around the active site was included. All QM/MM computations made use of density functional theory (DFT) for the QM part and the CHARMM force field for the MM part. They employed the ChemShell software38 using the TURBOMOLE program suite39 for QM, and the DLPOLY40 code for MM calculations. Geometry optimizations were carried out with the hybrid delocalized internal coordinate (HDLC) optimizer.41 In these calculations all the atoms except those contained in a sphere of r ) 10 Å around the inhibitor molecule were kept fixed. The convergence criterion was set to the default value (0.00045). Classical MD runs were performed to obtain additional information about the reversible intermolecular interactions which stabilize the reactants (enzyme-inhibitor complex before the ring-opening reaction) and the products (enzyme-inhibitor complex after the ring-opening reaction) of the irreversible inhibition step. These investigations started from the structures obtained in QM/MM geometry optimizations. The water sphere was again expanded to 50 Å around the enzyme, and the resulting structures were prepared using the same preparation procedure as described for the QM/MM. The subsequent MD run involved 15 ps heating, 200 ps equilibration, and a 200 ps sampling in an NVT system using the Berendsen thermostat.42 In all MD runs the Verlet leapfrog algorithm43 was used combined with a time step of 1 fs. In the case of the reactants, the interactions between the sulfur atom of the Cys29 and the attacked carbon atoms pose a problem for the MM treatment, since the bond-forming process has already started to some extent. This interaction was switched off in these classical MD simulations by setting the charge of the sulfur center to nearly zero (-0.07). If this is not done, conformations with hydrogen bonds between sulfur and the various hydrogen-bond donors of the inhibitor (e.g., N6H) are sampled. They may well be adopted in the initial stages of the inhibition process, but they are irrelevant to us since these geometrical arrangements cannot lead to the irreversible inhibition reaction studied presently. Also for these MD runs the CHARMM force field was used. To describe the epoxide some parameters concerning the ring structure had to be introduced. We used charges determined in the QM calculations and bond and angle parameters corresponding to sp2-hybridized C atoms. For the O atom of the ring the parameters of an ester oxygen were employed. Equilibrium distances and angles were set to the values obtained in QM computations. The QM part (Figure 1) generally contains the electrophilic warhead of the inhibitor with different substituents (see below)

Mladenovic et al. and the active site of the cysteine protease consisting of His199 and Cys29 residues (cathepsin B numbering).44 On the basis of many experimental results10 the zwitterionic form (Cys29-S-/ His199-H+) was adopted. In our standard ansatz the protein residues that form the oxyanion hole are not included in the QM part. Their interaction with the active site is mainly of electrostatic nature and should thus be well represented by the electrostatic embedding approach. To validate this assumption we performed additional computations in which the oxyanion hole was included in the QM part. Reaction energies obtained for this large QM region differed less than 0.3 kcal mol-1 (0.7%) from the results obtained with our standard ansatz, which proves its accuracy. To differentiate the influence of the warhead and of the peptidomimetic side chain, the inhibitor was first modeled by the three-membered ring together with the carboxylate substituent, and water molecules were inserted into the space occupied by the side chain (during the preparation phase). In the following this system will be denoted as OX-COO-. To study possible protonation processes accompanying the ring-opening reactions23,32 an ammonium molecule was included into the QM part as a proton source. The influence of the peptidomimetic side chain was studied by employing the actual E64c inhibitor (rather than water molecules) in the calculations. This system will be called E64c. In these computations only the warhead and the carboxylate substituent were included in the QM part, whereas the peptidomimetic side chain was treated at the MM level. In the E64c system the ammonium ion was omitted. Aziridine-based inhibitors were modeled by a system which is analogous to the OX-COO- system (AZR-COO-; AZRHCOO-, see below). To investigate the importance of the His199 inhibitor interaction on the overall reaction profile, the computations were performed both with a protonated or deprotonated ND1 center of the His199 residue (Figure 1); they will be distinguished by the terms protonated His199 and deprotonated His199. The investigated inhibition reactions involve the SN2 ring opening of the inhibitor (mainly described by angle ∠OCRCβ or ∠OCβCR) and the creation of the Scys29-C bond (mainly described by the distance RCR-S or RCβ-S). Due to their interplay both internal coordinates had to be used as reaction coordinates leading to a two-dimensional potential energy surface (PES). For each point of the PES corresponding to two fixed values of reaction coordinates (R-attack: ∠OCβCR and RCR-S; β-attack: ∠OCRCβ and RCβ-S) all other internal degrees of freedom were optimized to find the minimum energy path of the inhibition process. Points on the PES close to the transition state served as input for a TS search with the HDLC optimizer. Thereafter the optimized TS structures were relaxed to obtain the reactants and products to ensure that the stationary points are connected by a contiguous path. The geometries of the stationary points were used for the determination of barriers and reaction energies. During geometry optimizations and PES scans QM calculations were done at the RI45-DFT/B-LYP46/TZVP47 level of theory, whereas energies were determined by single-point computations at the DFT/B348-LYP/TZVP level. DFT is wellknown to describe many properties with an excellent costbenefit value.49-51 Nevertheless, this is often based on an error compensation which does not work in all cases,52-55 and indeed in many cases multireference approaches are necessary to obtain reliable reaction paths.56-59 For the present problem, however, our procedure was found to be sufficiently accurate for the determination of the PES.31 Problems encountered for protonated epoxides60 are not relevant presently since the protein environ-

QM/MM Calculations of Epoxide-Based Inhibitors

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5461

ment is not sufficiently acidic to protonate the reactants. Protonation only occurs far beyond the TS so that the kinetics are not influenced. The products are protonated alcohols for which similar problems are not known. In cases where the B-LYP-based PESs were too flat to determine suitable starting structures for the TS search, B3LYP single-point calculations at the RI-DFT/B-LYP structures were performed around the supposed TS area to obtain an appropriate starting point for the TS search at the DFT/B3-LYP/ TZVP level. This procedure was normally sufficient for locating the saddle point. Further optimizations to the product and reactant were performed at the same level of theory. 3. Results Epoxide-Based Inhibitors. To determine the underlying effects governing the observed inhibition potencies of epoxidebased inhibitors, the PESs of the ring-opening reactions with the thiolate moiety of the Cys29 residue were computed. Figure 1 sketches the studied reactions and denotes the important geometrical parameters and the enumeration of the atomic centers. The inhibitor was represented either by an epoxide ring with a carboxylate substituent (OX-COO-) or by the complete E64c inhibitor (E64c). To investigate the regioselectivity the attack of the thiolate group of the Cys29 residue at the C2 (Figure 1: R-attack) and at the C3 center (β-attack) were computed. The importance of the ionic interaction between the inhibitors’ carboxylate and the active site histidinium ion was assessed by calculations in which the ND1 center of the His199 residue (Figure 1) is protonated (abbreviated as protonated His199) or deprotonated (deprotonated His199). The computed PESs for the irreversible inhibition step are given in Figure 2 for the OX-COO- system and in Figure 3 for the E64c (in each case with protonated His199). They clearly show that both the opening angle ∠OCβCR (∠OCRCβ) and the distance between the sulfur and attacked carbon center, RCR-S (RCβ-S), contribute to the reaction coordinate. For the R-attack at the OX-COOsystem ∠OCβCR and RCR-S change similarly from the very beginning of the reaction. For the corresponding β-attack the change in ∠OCRCβ starts somewhat delayed. For the E64c both coordinates change from the very beginning for both reactions. The PESs obtained for the deprotonated His199 residue are qualitatively similar and are therefore not shown. Figure 4 summarizes the computed reaction barriers and energies. The products PH and P differ in the protonation of the alcoholate group which occurs through a proton transfer from the nearby NH4+ cation. As already noted before23,31 such protonation takes place far behind the transition state (TS) of the ring-opening reaction. The TS remains unprotonated for epoxide-based inhibitors even for proton donors stronger than NH4+. Selected geometrical parameters of the stationary points for the OX-COO- and the E64c model reactions are listed in Table 1 which contains both the hydrogen-bond distances and the corresponding heavy-atom distances, for comparison with X-ray data. The computed structures of the products of the reactions with E64c obtained with a protonated ND1 center of the His199 residue are shown in Figures 5 and 6. Our previous computations have predicted these reactions to be exothermic already for the nonsubstituted61 case. This is now also found for the R- and the β-attack at all substituted epoxides (Figure 4) indicating that the experimentally determined regioselectivity is a kinetic effect and that the irreversibility of the inhibition is not influenced by substituents. By contrast, substituent effects can substantially change the kinetics. Let us first concentrate on situations in which the ND1 center of the

Figure 2. PES obtained in QM/MM calculations for the OX-COO-R system (top) and the OX-COO--β system (bottom). The CR atom of the reactant has (S) configuration.

His199 residue is protonated (Figure 4, Table 1, entries 1-12). For the reaction of an unsubstituted epoxide ring with cathepsin B a reaction barrier of 17 kcal mol-1 has been predicted for the R-attack (Figure 4).61 This compares to a barrier of about 14 kcal mol-1 computed for the reaction of an unsubstituted epoxide with a free thiolate in a polar solvent. For such a solvent environment a carboxylate substituent raises the barrier of the R-attack to about 17 kcal mol-1, whereas the reaction for the β-attack remains at about 14 kcal mol-1.32 For the enzyme environment both barriers drop dramatically showing the strong influence of the enzyme-inhibitor interactions on the ringopening reaction. The resulting QM/MM barriers are only 2 kcal mol-1 for the R-attack and 9 kcal mol-1 for the β-attack (Figure 4, left-hand side: OX-COO-). The enzyme environment thus favors the R-attack and determines the regiospecificity found by X-ray measurements. The geometrical arrangements of the stationary points (Table 1, entries 1-3) reveal the origin of the strong influence of the enzyme environment on the kinetics. A very strong hydrogen

5462 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Figure 3. PES obtained in QM/MM calculations with the E64c inhibitor. Top, R-attack; bottom, β-attack.

bond is present between O18 and HD1 (Table 1: RO18-ND1 ) 1.6-1.7 Å). Additionally, slightly weaker hydrogen bonds are found between the carboxylate and the oxyanion hole, i.e., between O17 and the NH2 group of the Gln23 residue (RO17-NE2 ≈ 1.9 Å) and between O17 and the NH group of the Cys29 residue (RO17-NCys29 ) 2.0-2.3 Å). Judging from these distances the hydrogen-bond network is quite strong and should thus play an important role in the inhibition process. A comparison between the R- and the β-attack (Table 1, entries 1-3 vs 4-6) shows that RO18-ND1 and RO17-NE2 are nearly identical, whereas RO17-NCys29 elongates during β-attack, implying that the necessary motion of the sulfur atom to the more distant C3 atoms weakens this interaction to some extent. Another reason of the very low barrier for the R-attack becomes clear when the geometrical parameters of the reactants are compared. For the reaction of an unsubstituted epoxide with the active site of cathepsin B previous studies61 gave RC2-S ≈ 4 Å and -OC3C2 ) 60° for the reactants but RC2-S ) 2.8 Å and -OC3C2 ) 67° for the TS. For the R-attack at a carboxylated epoxide the present study yields RC2-S ≈ 3 Å and

Mladenovic et al. -OC3C2 ≈ 65° for the reactants and RC2-S ≈ 2.8 Å and -OC3C2 ) 71° for the TS (Table 1, entries 1, 2), i.e., the reactant geometries are strongly influenced, but the transition structures are quite similar. This indicates that the strong interactions between the carboxylate moiety and the enzyme pull the inhibitor into the active site toward the attacking thiolate group of the Cys29 residue such that the ring already starts to open in the reactant valley. In other words, the strong interactions within the hydrogen network push the reactants along the reaction path uphill toward the TSs leading to considerably smaller barriers. This effect also influences the regiospecificity, i.e., the difference in the reaction barriers of the R- and the β-attack. In the corresponding reactants, the R-attack angle -OC3C2 has a value of 65°, whereas the β-attack angle -OC2C3 is 59°, i.e., even slightly smaller than the value of 60° found for a solvent environment. The interaction between the thiolate of the Cys29 residue and the three-membered ring that is already present in the reactants thus prepares the ring more for the R- rather than the β-reaction. To study this in more detail, Figure 7 plots the correlation between the computed reaction barriers and the corresponding variation in the ring-opening angle ∆-OCC ) -OCC (TS) - -OCC (reactant). Obviously there is such a connection for all epoxide- and aziridine-based inhibitors studied presently. A similar correlation is also found for the relative change in the RC2-S distance (Figure 7: bottom). This influence on the regioselectivity can be summarized as follows: When passing from the reactant to the transition structure, the reaction coordinate of the R-attack, -OC3C2, has to increase by only 6°, whereas the corresponding angle for the β-attack (-OC2C3) has to change by 17°, from 59° to 76°. This larger distortion will raise the barrier for the β-attack which is computed to be about 7 kcal mol-1 higher than that for the R-attack. The peptidomimetic side chain has only a minor influence on the kinetics of the R-attack. There is only a small difference in the activation barriers computed for OX-COO- and E64c, for which a value of 1 kcal mol-1 is predicted (Figure 4, righthand side). It should be noted, however, that the interactions between the side chain and the enzyme influence the hydrogenbonding network connected with the carboxylate. Again the O17-NCys29 hydrogen bond is affected most, since its length in the reactants and TSs increases by about 0.5-0.6 Å when going from OX-COO- to E64c. Similar but considerably smaller variations are also found for RO17-NE2, whereas RO17-NE2 remains essentially unchanged. These two other hydrogen bonds ensure that the inhibitor is still pushed toward the attacking thiolate moiety which explains the small overall effect on the barrier height of the R-attack. Concerning the interactions of the peptidomimetic side chain itself, the O5-Gly74 hydrogen bond seems to be the most important. The distance between N6H and the Gly198 residue is already considerably longer (3.8-4.2 Å), and the heavy-atom distances within the hydrophobic pockets are even longer (about 5 Å), thus implying rather weak interactions. Due to compensation effects the peptidomimetic side chain is thus less important for the kinetics of the R-attack, but it considerably influences the regioselectivity. The computed barriers for the β-attack (attack at C3, Figures 1 and 4) are 7 (OX-COO-) and 15 kcal mol-1 (E64c) higher than the corresponding barriers of the R-attack. This shows that the 100% regiospecificity found in X-ray measurements results equally from the hydrogen network connected with the carboxylate substituent on the one hand and from interactions of the

QM/MM Calculations of Epoxide-Based Inhibitors

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5463

Figure 4. Reaction profiles computed for the reactions of epoxide-based inhibitors with respect to the reactants R. For definition of the various simulations see text. TS denotes the transition states, and P and PH denote the deprotonated and protonated products, respectively.

TABLE 1: Selected Geometrical Parameters for the Stationary Points of the Reaction Profiles Computed for the Reactions of Epoxide-Based Inhibitorsa entry

Rx

∠δ

RO18-ND1

RO17-NE2

RO17-NCys29

RO5-Gly74

RN6-Gly198

-

-

1 2 3

R TS PH

3.07 2.81 1.92

64.7 71.1 112.2

OX-ΧΟΟ-/R-Attack/Protonated His199 1.66/2.72 1.94/2.91 2.31/3.25 1.62/2.62 1.92/2.88 2.20/3.14 1.58/2.65 1.88/2.84 1.99/2.94

4 5 6

R TS PH

3.33 2.76 1.88

59.2 76.1 108.3

OX-COO-/β-Attack/Protonated His199 1.67/2.72 1.95/2.92 2.51/3.42 1.66/2.72 1.90/2.87 2.52/3.42 1.59/2.66 1.87/2.84 2.21/3.12

7 8 9

R TS P

3.07 2.80 1.91

64.0 72.4 111.0

E64c/R-Attack/Protonated His199 1.65/2.68 2.21/3.15 1.64/2.65 2.15/3.10 1.67/2.64 1.91/2.88

2.97/3.85 2.67/3.59 1.96/2.91

2.00/2.94 1.99/2.93 2.01/2.95

3.88/4.02 3.69/4.03 3.75/4.05

10 11 12

R TS P

3.56 3.06 2.01

57.2 79.6 107.4

E64c/β-Attack/Protonated His199 1.66/2.70 2.38/3.29 1.66/2.70 2.28/3.21 1.62/2.63 1.96/2.93

3.40/4.25 3.30/4.15 2.64/3.54

2.10/3.03 2.08/3.01 2.00/2.93

3.51/3.81 3.49/3.82 3.58/4.06

13 14 15

R TS P

3.41 3.14 1.92

63.5 70.7 111.1

E64c/R-Attack/Deprotonated His199 -/4.38 2.26/3.14 -/4.30 2.28/3.15 -/4.20 2.17/3.05

2.13/2.93 2.08/2.91 1.80/2.73

1.91/2.85 1.91/2.85 2.01/2.93

2.91/3.83 3.18/4.00 3.73/4.28

a Rx gives the distance between the sulfur center of the thiolate and the attacked carbon atom. ∠δ denotes the angle describing the ring opening (O1C3C2 for the R-attack and O1C2C3 for the β-attack). All distances are given in angstroms, angles in degrees. For definition of the various parameters see Figures 1 and 6. The last five columns refer to hydrogen-bond geometries (first value, X‚‚‚H distance; second value, heavy-atom distance).

peptidomimetic side chain with the S1-S3 subsites on the other hand (Figure 6). How both types of interactions affect each other, and thereby the regioselectivity, can be best seen from the hydrogen-bond network connected with the carboxylate: the O17-NCys29 bond shows again the strongest differences when R- and β-attack are compared (Figure 6, Table 1). For the OX-COO- system we have already pointed out that the interaction between the thiolate group of the Cys29 residue and the C2 center prepares the three-membered ring for an Rrather than a β-reaction. This effect is even more pronounced in the E64c system where the ∆-OCC values are 22.4° for the β-attack and only 8.4° for the R-attack. The increased difference in comparison to the OX-COO- system (6.4° vs 16.9°) nicely correlates with the increased difference in the barrier height. This shows that all effects leading to the regioselectivity (direct reversible interaction between inhibitor and enzyme, preference of the R-reaction due to interactions between the negatively charged sulfur and C2) are enhanced when going from the OXCOO- to the E64c system.

The small barrier computed for the irreversible step in the enzyme questions the commonly accepted two-step model for such inhibitors. Due to the small reaction barrier of the ringopening reaction the irreversible step can be expected to take place immediately for suitable geometrical arrangements. Consequently, the assumed equilibrium of the reversible enzymeinhibitor complex and the separated compounds will be mainly determined from conformations of the inhibitor which do not allow the ring-opening reaction, e.g., with the thiolate too far away from the carbon centers of the ring. According to our calculations, for suitably oriented conformations the two reaction steps will effectively merge into a single one, and hence the first-order rate constant of inhibition, ki, is not only affected by the irreversible ring-opening reaction but also by the reversible network interactions. In the commonly accepted two-step model the latter interactions are subsumed under the dissociation constant of the noncovalent enzyme inhibitor complex, Ki. In our scenario an absolute discrimination between ki and Ki is not possible. For example, if Ki is increased due to more

5464 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Mladenovic et al.

Figure 5. R- (top) and β- (bottom) product structures of the E64c inhibitor. The residues of the active site and the most important interactions are shown. Dots indicate the hydrophobic pocket.

stabilizing interactions in the hydrogen-bonding network, ki, concomitantly increases. Nevertheless, it is still helpful to distinguish between effects which influence the formation of the reversible enzyme-inhibitor complex and interactions which affect the ring opening. Our results indicate that the secondorder rate constant of inhibition, k2nd, which includes the effects resulting from the reversible and the irreversible interactions, is more suitable for characterizing this class of irreversible inhibitors. QM/MM geometry optimizations provide a static picture. To check if the flexibility of the system influences the importance of specific interactions, classical MD simulations for reactants

and products were performed. The geometries predicted by the QM/MM procedure were taken as starting points. Figure 8 shows the distances between the sulfur center of the Cys29 thiolate and the attacked carbon center (black: CR, red: Cβ) along the MD trajectory. They vary between 3.3 and 4.2 Å (average 3.6 Å) for the R-attack and between 3.5 and 4.5 Å (average 4.0 Å) for the β-attack, consistent with the preference of the R-attack. The fluctuations of the (hydrogen-bond) distances along the trajectories of classical MD simulations allow an assessment of the strength and importance of the underlying interactions. The MD results for selected coordinates are plotted in the

QM/MM Calculations of Epoxide-Based Inhibitors

Figure 6. Sketch of the network which stabilizes the inhibitor by reversible interactions. The product of the R-attack is given.

Supporting Information. They nicely agree with the conclusions derived from the optimized QM/MM distances. In the reactants, the interactions between the carboxylate moiety of the inhibitor and His199 (averaged distances RO18-ND1 (av) ≈ 2.7 Å, variations