CHEMPHYSCHEM ARTICLES DOI: 10.1002/cphc.201402542

Protocol for Rational Design of Covalently Interacting Inhibitors Thomas C. Schmidt,[a] Armin Welker,[b] Max Rieger,[c] Prabhat K. Sahu,[d] Christoph A. Sotriffer,[b] Tanja Schirmeister,[e] and Bernd Engels*[a] The inhibition potencies of covalent inhibitors mainly result from the formation of a covalent bond to the enzyme during the inhibition mechanism. This class of inhibitors has essentially been ignored in previous target-directed drug discovery projects because of concerns about possible side effects. However, their advantages, such as higher binding energies and longer drug-target residence times moved them into the focus of recent investigations. While the rational design of non-covalent inhibitors became standard the corresponding design of covalent inhibitors is still in its early stages. Potent covalent inhibitors can be retrieved from large compound libraries by covalent docking approaches but protocols are missing that can

reliably predict the influence of variations in the substitution pattern on the affinity and/or reactivity of a given covalent inhibitor. Hence, the wanted property profile can only be obtained from trial-and-error proceedings. This paper presents an appropriate protocol which is able to predict improved covalent inhibitors. It uses hybrid approaches, which mix quantum mechanical (QM) and molecular mechanical (MM) methods to predict variations in the reactivity of the inhibitor. They are also used to compute the required information about the noncovalent enzyme–inhibitor complex. Docking tools are employed to improve the inhibitor with respect to the non-covalent interactions formed in the binding site.

1. Introduction The residence time of a given drug bound to its target is one of the key figures for its efficiency. It depends on many factors but the interaction strengths between inhibitor and enzyme is one of the essential parameters. In that respect non-covalently interacting inhibitors possess disadvantages because their binding free energies to enzymes seldom exceed 60 kJ mol 1.[1 2] Covalent inhibitors can achieve higher binding energies, but based on early works on bromobenzene and acetaminophene[3–7] the general chemical reactivity of some classes of covalent inhibitors has led to concerns against covalent inhibitors in general.[8] Hence, despite various very successful examples of covalent drugs (e.g. penicillins, acetylsalicylic acid) covalent

inhibitors are avoided when starting a target-directed drug discovery project. It is thus not surprising that essentially all existing first-in-class covalent drugs result from biological assays and that their covalent molecular mechanism was found later on. The possibility to achieve higher binding energies has led to a certain resurgence of covalent inhibitors[1, 8] but shortcomings in their rational design hamper their development. For the discussion of the corresponding approaches the inhibition mechanism of covalent inhibitors has to be considered in some detail. This is depicted in Figure 1 together with a sketch of the corresponding energy diagram. In the first step, enzyme E and inhibitor I form an enzyme–inhibitor complex whose binding free energy DGb only depends on non-covalent interac-

[a] T. C. Schmidt, Prof. Dr. B. Engels Institut fr Physikalische und Theoretische Chemie Julius-Maximilians-Universitt Wrzburg, Wrzburg (Germany) E-mail: [email protected] [b] A. Welker, Prof. Dr. C. A. Sotriffer Institut fr Pharmazie und Lebensmittelchemie Julius-Maximilians-Universitt Wrzburg, Wrzburg (Germany) [c] M. Rieger Max Planck Institut Karlsruhe (Germany) [d] P. K. Sahu Department of Chemistry National Institute of Science and Technology Palur Hills, Berhampur-761008 (India) [e] Prof. Dr. T. Schirmeister Institut fr Pharmazie und Biochemie: Therapeutische Lebenswissenschaften Johannes Gutenberg-Universitt Mainz (Germany) Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/cphc.201402542.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Figure 1. Two-step inhibition mechanism of covalent inhibitors. Step1: formation of non-covalent enzyme–inhibitor complex E···I. Step 2: chemical reaction leading to the covalent enzyme–inhibitor complex E–I.

ChemPhysChem 2014, 15, 3226 – 3235

3226

CHEMPHYSCHEM ARTICLES tions (E···I). In the second step, the chemical reaction that leads to the covalent enzyme–inhibitor complex E–I takes place. The accompanied additional contribution to the binding free energy of the enzyme inhibitor complex is the free reaction energy of the chemical reaction (DGR). If it is so exergonic that no back-reaction can take place, the enzyme is irreversibly blocked (covalent irreversible inhibitor). If the reaction is only weakly exergonic the inhibition is covalent reversible because the back reaction is still possible. For endergonic reactions no or nearly no additional covalent binding effects are found. Most molecular docking procedures focus on the docking of non-covalent inhibitors. Possibilities to retrieve covalent inhibitors from compound libraries are offered for example by the packages GOLD,[9, 10] FlexX[11] or QXP/FLO.[12] Such approaches are named “Covalent Docking”. Such implementations have been used with some success, however, several drawbacks seem to exist.[13, 14] Selzer and co-worker used the implementation in the GOLD package to retrieve potent covalent cathepsin K inhibitors in a high-throughput range from large compound libraries.[13] The study focused on compounds possessing an electrophilic carbon centre which reacts with the sulphur active site (Cys 25). Selzer and co-worker indeed found three inhibitors with nanomolar IC50 values, but the success seems to result mainly from a manual selection step. The original automated docking ranks of the three compounds (out of 1247 compounds) were 242, 334, and 621 if GOLDScore was used as scoring function. The scoring function ChemScore predicted the rankings 322, 571, and 616. Using cognate covalent ligand docking to test the accuracy of this approach, Selzer and co-worker noticed that “the scoring functions in general, are prone to errors in which the poses closest to the crystal structure conformation do not consistently receive the best score.”[13] In the algorithm implemented in GOLD the ligand is docked in the active site under the restriction that the link atom between ligand and enzyme keeps the position of the sulphur centre of the cysteine unit. To take those geometrical changes in the ligand into account that arise from the formation of a covalent bond to the enzyme, a methylthiolate moiety was added to the electrophilic centre of the ligand. Subsequently, the resulting unit was geometry-optimized employing the MMFF94 force field. It is important to note that this approach is applied to the covalent enzyme–inhibitor complex, that is, the final stage of the inhibition mechanism (Figure 1). The implementation in Autodock uses Gaussian wells in the scoring function to take the covalent bonding into account. Problems with this approach seem to arise because these Gaussian wells have to be set to specific reaction coordinates and because the necessary parameters had to be carefully optimized to achieve reliable predictions.[14, 15] To increase the accuracy Kwoh and co-worker replaced the Gaussian wells by Morse potentials, which were carefully parameterized on the basis of ab initio results.[14] Using this approach they achieved considerably higher success rates in cognate covalent docking than the original Autodock approach. Please note that also these approaches simulate the stability of the covalent enzyme–inhibitor complex.  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org Some errors in the predictions of covalent docking approaches described above may result because they focus on the covalent complex E–I, while effects connected to the noncovalent complex E···I (formed in the first step) are largely neglected. However, the stability and geometry of E···I strongly influence the efficiency of the subsequent chemical reaction. If the non-covalent complex has a very short lifetime only very fast reactions can take place. Slower reactions with higher reaction barriers cannot take place because the complex is already decomposed before the chemical reaction has started. The geometry of the non-covalent complex is important because the subsequent reaction can only happen if the reaction centres of ligand and enzyme are sufficiently close to each other. The strong influence of this effect was shown to be the reason for different inhibition potencies of the various stereoisomers of E64c for cathepsin B.[16] It is also the underlying basis for the regioselectivity of this compound in the reaction with cathepsin B.[16] Both studies indicate that already small increases in the average distance between the reaction centres diminish the reaction probability even for strong exothermic reactions. Taking this into account, retrieval procedures of potent covalent inhibitors should also consider the properties of this complex. The consideration of the non-covalent complex is even more important for the rational design of improved covalent inhibitors. It should start with a step in which the ligand is varied in order to improve the binding free energies for the non-covalent inhibitor enzyme complex. At the same time the orientation of the ligand in this complex should be optimized with respect to the targeted inhibition reaction. It is of course important to consider the influence of these variations on the subsequent covalent bond formation step because the reaction might be hampered by the modifications introduced. However, this approach contains some difficulties. In general, no data exist for the non-covalent complex. Most available crystal structures represent the situation in the covalent inhibitor–enzyme complex because the lifetime of the non-covalent counterpart is too short due to the inherent reactivity of the warheads. The second problem arises because a detailed computation of the second step deserves quantum chemical approaches since it involves bond forming and breaking processes. This paper presents an appropriate protocol for the rational design of covalently interacting inhibitors (Figure 2). It uses QM/MM hybrid approaches, which combine quantum mechanical (QM) and molecular mechanical (MM) methods, for the description of the chemical reaction (Figure 1: step 2). Such methods are also used to extract information about the non-covalent complex E···I from the crystal Figure 2. Protocol for develstructure of the covalent complex opment and improvement of E–I. Modifications of the inhibitor covalent inhibitors. ChemPhysChem 2014, 15, 3226 – 3235

3227

CHEMPHYSCHEM ARTICLES are then achieved by applying standard docking tools to the non-covalent complex E···I. Their influences on the chemical reactivity of the inhibitor are subsequently tested by further QM/ MM computations. The protocol is developed to improve the inhibition potency of a given inhibitor. Our approach is an addition to the covalent docking approaches which focus on the retrieval of potent covalent inhibitors from large compound libraries. Our approach can be used subsequently to further improve promising inhibitors detected by such approaches. To describe the capabilities of our approach (see Figure 2), we have chosen HIV-1 PR as test case, because an overwhelming number of inhibitors is known and corresponding crystal structure datasets are available. We used the crystal structure of the covalent complex with 1,2-epoxy-3(p-nitrophenoxy)propane (EPNP, Figure 3) as starting point for our investigation

Figure 3. Chemical structures of EPNP and AZNP.

and varied recognition unit and warhead to predict improved covalent inhibitors. EPNP is a non-peptidic inhibitor of pepsin[18] and other aspartic proteases such as HIV-1 PR.[19] It has been shown to react in a stoichiometry of one molecule per dimeric enzyme, suggesting the esterification of either of the two catalytically active Asp 25 residues.[20] The reaction scheme of the underlying bond formation is depicted in Figure 4. The inhibition constant Ki was determined to 9.9  1.0 mm to 11  2 mm with an inactivation rate of 0.04  0.005 min 1 to 0.06  0.006 min 1, depending on the reference.[19, 20]

1.1. Protocol for Rational Design of Covalent Inhibitors The cycle developed for the rational design of covalently interacting inhibitors is sketched in Figure 2. Our protocol starts

www.chemphyschem.org with the selection of an appropriate crystal structure of a covalent enzyme–inhibitor complex. It is helpful if the chosen crystal structure already largely resembles the targeted inhibitor, but this is not a prerequisite. The required structural information about the non-covalent enzyme–inhibitor complex E···I is obtained in a preparatory step of our protocol. To gain such information, we use QM/MM approaches (see below) to compute the reaction path backwards starting with E–I for which a crystal structure is available and ending at E···I. QM/MM hybrid approaches are necessary because quantum mechanical approaches (QM) are needed to describe the formation and breaking of covalent bonds while molecular mechanical (MM) methods are needed to take the influence of the enzyme environment into account. Please note that pure QM approaches are far too expensive due to the size of the involved enzymes. After computing the necessary structural information of the non-covalent enzyme–inhibitor complex the cycle of the rational design can be started. Variations of the inhibitor to improve its affinity towards the target are obtained by applying standard docking tools to E···I. They also allow an estimate to what extent the stability of the complex increases. It is of course also possible to vary the warhead or its substitution pattern. Such modifications may not only influence the affinity, but also change thermodynamics and kinetics of the second step of the inhibition mechanism. The structural information obtained for the modified E···I complex in this step is then refined by MD simulations to obtain a starting structure for subsequent QM/MM computations. In these subsequent QM/MM computations we test how the modifications of the inhibitor influence thermodynamics and kinetics of the second step of the inhibition mechanism. If they are favourable, the cycle can be repeated to obtain further improvements of the new inhibitor model. If they are unfavourable, we have to go back again to develop different variations of the parent compound again using docking. Our protocol should provide quite reliable predictions due to three reasons.[16, 17, 21–24] The geometries of E–I with which our rational design protocol starts are quite accurate because they are taken from crystal structure data. The same holds true for the geometries of E···I which are computed via QM/MM. Our confidence results from the observation that QM/MMbased refinements of crystal structures were already proven to provide more accurate results than the standard approaches.[23, 25–27] The computed kinetics and thermodynamics

Figure 4. Reaction scheme of the bond formation between EPNP and the active site Asp 25 of HIV-1 PR.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

ChemPhysChem 2014, 15, 3226 – 3235

3228

CHEMPHYSCHEM ARTICLES of the reactions (starting and modified inhibitor) will also be reliable.[22] Using comparable approaches as used here, we were, for example, able to clarify the inhibition mechanism of cathepsin B by E64c.[23, 25] These calculations could explain the regio- and stereoselectivity of the inhibition.[16, 17] Furthermore, suggestions for improved inhibitors based on QM/MM calculations were proven experimentally.[28]

Computational Methods Classical molecular dynamics simulations (MD) were performed using the NAMD program package[29] employing the CHARMM force field.[30] The protein was modelled using the CHARMM 27.1 force field, water molecules were treated according to the TIP3P parameterization.[31] Inhibitor molecules were parameterized via the web based interface of SWISSPARAM[32] and the corresponding parameters were refined roughly following the CHARMM General Force Field[33] procedure. For simulations on inhibitors covalently bound to the active site Asp 25, CHARMM patches have been generated linking the OD2 atom of the aspartate to the electrophilic carbon atom of the inhibitor while cleaving the adjacent bond to the heteroatom and adding a hydrogen atom at the liberated alcoholate/amide. To model aqueous solution, the systems were solvated by a water shell of 60  radius centered at the protein’s centroid. The non-covalently complexed inhibitors of the (1-(1,3-diphenylpropyl)-2-ethyl-3-{[(4-nitronaphthalen-1-yl)oxy]methyl}aziridine) group shown in Figure 5 were modelled as neutral molecules with

www.chemphyschem.org For the QM/MM calculations, the ChemShell[36–38] package was used employing the internal DL POLY force field to compute the molecular mechanical, and Turbomole[39] for the quantum mechanical subsystem. The interactions are represented within the electrostatic embedding Scheme while sliced bonds are treated by link atoms to avoid issues due to incomplete valence shell populations at the QM/MM boundary. The QM part is calculated on density functional theory (DFT) level employing the BLYP[40, 41] functional for the potential energy profiles. All QM calculations are performed using the TZVP[42] basis set available with Turbomole. To reduce the computational costs during the QM/MM calculations, the original water shell is truncated at a radius of 15  around the active site residues while the protein is still modelled completely. All residues outside a 10  radius are fixed with the remaining sphere being the active region of the QM/MM computations. The QM part comprises up to 50 atoms from both catalytic aspartate residues, truncated between the a-carbon atom and side chain, the epoxide or aziridine ring or 3,6-dimethyl-3,6-diazabicyclo[3.1.0]hexane, respectively, for the inhibitor. Water molecules were also included if they were needed to model mediated proton transfers. Docking experiments using the GOLD SUITE 4.3 from the CCDC with standard configuration[9, 43] were performed to modify the affinity of the ligand to the binding pocket of HIV-1 PR by modification of the substitution pattern of the functional group. The employed scoring functions were GoldScore, ChemPLP and SFCscore.[44] In case of SFCscore, the sfc 290m function was used. The presented scores are the ones with water set up as cofactor. Without water, the differences were usually smaller than 0.1 pKi. Poses for modified inhibitor models were rescored using ChemPLP, GoldScore and SFCscore.[44]

2. Results and Discussion 2.1. Step1: Preparation of the Starting Structures for Docking

Figure 5. Chemical structure of the DPPA inhibitors derived from the Michaelis complex of AZNP.

the active-site dyad residing in the expected protonation state of the free enzyme. In the case of inhibitors based on the 3,6diazabicyclo[3.1.0]hexane scaffold (1020 Figure 6), several protonation states have been considered. Depending on the environment, either the pyrrolidine (pKa of the corresponding protonated species  11.27)[34] , the aziridine nitrogen centre (pKa  7.9)[35] or both are expected to be protonated. MD simulations were carried out at a temperature of 310 K up to a simulation time of 5 ns. All simulations were preceded by an equilibration of the randomly assigned, Boltzmann-distributed velocities by means of a spherical harmonic potential, constraining each atom to its initial position. The corresponding force constant was halved every 5 000 steps and the equilibration was considered complete after 6 cycles.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Our protocol (Figure 2) starts with the selection of an appropriate crystal structure of the covalently bonded enzyme–inhibitor complex E–I. Due to the lack of a crystal structure containing HIV-1 protease complexed with EPNP, the necessary starting structure was obtained through the combination of the two crystal structures 1YTH[45] and 2SAM.[46] The crystal structure 1YTH[45] comprises the free dimeric HIV-1 protease and crystal water but the EPNP inhibitor is missing. The crystal structure 2SAM[46] provides the structure of the SIV protease covalently connected with EPNP. As both enzymes have nearly identical active sites and also coincide to a large extent (see the Supporting Information, Figure S1), the combination of both structures should provide a reasonable model for the HIV-1 enzyme covalently attached to EPNP. To test the accuracy of our force field approach we first modelled the free enzyme by molecular dynamics (MD) simulations starting from the crystal structure 1YTH.[45] It is generally believed that the two Asp moieties behave like a diprotonic acid with one Asp being protonated (carboxylic acid) and the other deprotonated (carboxylate). A comparison of more than 40 crystal structures of HIV-1 Protease (wildtype and mutants)[45–65] in uncomplexed as well as complexed form revealed ChemPhysChem 2014, 15, 3226 – 3235

3229

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

pends on their relative protonation state. When treated in the expected monoprotonated state, the orientation of Asp 25 and Asp 25’ agrees nicely with the orientation given in the crystal structure 1YTH (distance Cg(Asp 25)-Cg(Asp 25’)  5.0 ). Visual inspection even revealed the frequent occurrence of a cyclic structure containing 11 atoms comprising the carboxylate and the carboxylic acid moiety which are bridged by a water molecule.[67] If both Asp residues are deprotonated, the repulsion due the negative charges pushes the carboxylate groups 1.0  [distance Cg(Asp 25) Cg(Asp 25’) up to 6.1 ] further away from each other. The opposite effect is observed if both residues are modelled as carboxylic acid moieties. In this protonation state the Asp residues form a more stable dimer due to the two hydrogen bonds between the functional groups [distance Cg(Asp 25) Cg(Asp 25’) down to 3.6 ]. Our simulations are in agreement with previous simulations.[66–68] In the next step we performed computations to achieve information about the non-covalent enzyme–inhibitor complex between EPNP and the HIV-1 protease which is the starting point for the subsequent docking investigations. Our efforts started with MD simulations for the coFigure 6. Chemical structures of the four best ranked docking results for 3,6-diazabicyclo[3.1.0]hexane based inhibitors. Fitness scores determined with Goldscore are given in parenthesis. valent enzyme–inhibitor complex E–I (Figure 7 a) of HIV-1 PR and EPNP, which were built up from the overlay of the crystal almost identical geometries of the protein backbone for all the structures 1YTH[45] and 2SAM.[46] It reflects the situation after structures (see Figure S1 and Table S1). Larger structural deviastep 2 of the inhibition mechanism (Figure 1), that is, the tions are only found in the flap regions, where the protein three-member ring of the epoxide moiety is opened and its adapts to accommodate the complexed molecules.[66–68] A oxygen centre is already protonated. In accordance with the visual inspection of the results obtained from our MD simulacrystal structures our simulations indicate that the complexations of the uncomplexed enzyme showed good agreement tion by EPNP does not strongly influence the conformation of with the crystal structures mentioned above. The RMSD values the binding site residues of the HIV-1 protease. Our simulations for the positions of the heavy atoms of the protein backbone also indicate that the position of the ring-opened EPNP is stacompared to 1YTH remained well below 1.5 . Thus, the simubilized by hydrogen bonds to the flap water. lations prove the general usability of CHARMM in describing The computed potential energy surface is depicted in Figthe protein structure. The resulting structure is depicted in Figure S3. The geometrical orientations of the involved residues ure S2. for the corresponding non-covalent complex are sketched in Our MD simulations indicate that the relative orientation of Figure 7 lower half. The computed potential energy surface the two catalytic residues, Asp 25 and Asp 25’, strongly de 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

ChemPhysChem 2014, 15, 3226 – 3235

3230

CHEMPHYSCHEM ARTICLES

Figure 7. Sketches of EPNP within the active site of HIV-1 PR.

predicts that the inhibition reaction of the epoxide warhead starts with an attack of the OD2 of Asp 25 at the C2 of EPNP. The reaction barrier is computed to be about 30 kJ mol 1. According to our calculations, the reaction energy is about 120 kJ mol 1. Because the active site is water accessible, one  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org can assume that the emerging carboxylate group of Asp’ is complexed by a water molecule. Taken this into account, the reaction energy even drops to about 150 kJ mol 1. This indicates that the reaction is very exothermic, which is in line with the irreversibility of the inhibition of the HIV-1 PR by EPNP. Besides this, our computations also indicated a second possible mechanism mediated by water. In this reaction the proton transfer from Asp’ to the emerging alcoholate moiety of EPNP is mediated by one water molecule. For this reaction a barrier of about 45 kJ mol 1 is predicted while the reaction energy is 155 kJ mol 1. The higher barrier of the water mediated proton transfer indicates that the proton is directly transferred. Nevertheless, both processes are easily possible at ambient temperatures. Having all relevant geometrical information about the noncovalent EPNP HIV-1 protease complex at hand, docking investigations with the goal to improve the affinity of EPNP could start. However, we decided to use its aziridine analogue, AZNP (Figure 3), for the docking because its nitrogen centre offers additional possibilities. Hence, we replaced the epoxide warhead by an aziridine warhead and computed the corresponding inhibition reaction to ensure that the reactivity of the aziridine warhead is also suited for the inhibition. This is indeed predicted by our computations. The corresponding potential energy surface is given in Figure S4 while the variations in the geometries can be taken from Figures S5, S6 and S7. For the reaction with a direct proton transfer we computed a reaction barrier of 20 kJ mol 1 while about 50 kJ mol 1 are predicted for the water mediated mechanism. The reaction energy is computed to about 130 to 150 kJ mol 1, depending on the complexation of the resulting carboxylate group of Asp’ by water molecules. ChemPhysChem 2014, 15, 3226 – 3235

3231

CHEMPHYSCHEM ARTICLES 2.2. Step 2: Docking Investigations to Achieve Inhibitors with Improved Affinities in the Non-Covalent Complex After computing all relevant geometrical information about the non-covalent AZNP HIV-1 PR complex we performed docking investigations (see Computational Details) to achieve inhibitors with improved affinities of AZNP with respect to the active site of HIV-1 PR. In a first attempt to improve the affinity of AZNP, we added hydrophobic residues. The nitrophenyl substituent was extended to a naphthalene system. The aziridine hydrogen was replaced by a flexible alkyl scaffold with two phenyl rings. This ligand ((2S,3R)-1-[(R)-1,3-diphenylpropyl]-2-ethyl-3-{[(4-nitronaphthalen-1-yl)oxy]methyl}aziridine) (DPPA2, Figure 5) has a similar binding orientation compared to EPNP while additional pockets are addressed. A binding mode percentage of over 50 % was achieved and the pose was scored with an pKi value of 6.70 by SFCscore. A similar inhibitor with only a methylene linker between the upper right phenyl ring and the aziridine (DPPA1 Figure 5) was scored with 6.79. The difference is small but this extra flexibility should relieve the ring opening during the covalent reaction. Compared to the SFCscore-estimated pKi values of AZNP (3.66) these affinities are already increased by a factor of 1 000. To investigate how the introduced variations in the structure of the inhibitor influence its orientation in the active site we performed MD simulations. As starting point for the MD simulations we used the structures predicted by the docking procedure (see the Supporting Information for Cartesian coordinates) and monitored the distance between the reacting centres, namely the oxygen centres OD1 or OD2 of Asp’ and the attacked carbon centre C2 (Figure 8) during the MD run. If this

www.chemphyschem.org orientation does not switch back, that is, the new formation seems to be lower in energy than the starting orientation. The relaxed reorientation prevents the reaction because the nucleophilic carboxylate moiety of Asp’ cannot reach its electrophilic counterpart (C2-centre of the aziridine ring) anymore. Hence, we expect that the new inhibitor shows lower inhibition potency because the formation of the covalent enzyme–inhibitor complex is hindered. We avoid further QM/MM computations to investigate the energetics of the ring-opening reactions of DPPA1 and DPPA2 because they are quite expensive and the outcome of the MD simulations seems to be unambiguous. Our findings that higher affinities can be accompanied with a suppression of the reaction, underline the importance of the last step of our protocol. Hence, we went back and constructed a new set of structures again using docking. To avoid a rotation of the aziridine ring, the ether function was replaced by a pyrrolidine ring leading to the bicyclic scaffold shown in Figure 6. To avoid issues due to steric hindrance, the ethyl substituent was shortened to either methyl or even a single hydrogen atom. The substituent at the aziridine nitrogen was only marginally modified. We introduced a carboxylic acid function, which is directed into the solvent and further stabilizes the orientation within the corresponding pocket. Additionally, the length of the linker adjacent to the nitrogen atom was varied. On the basis of a comparison with commercially available drugs, we replaced the nitro function to the side chain known from tipranavir.[69, 70] The additional changes led to even better scoring results. For the smallest and best-ranked model compound 1020 (Figure 6) the pKi value was estimated to 8.45, that is, the affinity of compound 1020 is 50 000 fold higher than the affinity of AZNP. The three other compounds of the actual set perform comparably. The predicted orientation of 1020 within the noncovalent E···I complex is sketched in Figure 9. 2.3. Step 3: QM/MM Computations to Characterize the Reactions of the Improved Inhibitor

Figure 8. Distances between the nucleophilic oxygen centres of Asp 25 and the C2 atom of DPPA during an MD simulation started from the best ranked docking pose.

distance becomes large the reaction is hampered. It shows that the distance first remains at values of around 4 , but after about 300 ps it jumps to more than 6 . This jump results because the nitrogen centre of the inhibitor orients towards the aspartic acid residues. The trajectory also indicates that the  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

MD simulations of the 1020 HIV-1 PR complex revealed that the distance between the reactive centres, OD2 and C2 (Figure 9) stays sufficiently small. Hence, we performed QM/ MM computations to characterize the reaction of 1020 with HIV-1 PR. In line with our results for the free enzyme, we assume that only the Asp 25’ residue is protonated. Under physiological conditions (pH  7.5) the nitrogen centre of the pyrrolidine unit (Figure 9 N1: pKa = 11.3) is also protonated. For the aziridine nitrogen centre (Figure 9 N2: pKa = 7.9) both possibilities have to be considered. To take into account that the active site is water accessible, additional water molecules were added. Their positions were taken from MD simulations. Figure 10 displays the two-dimensional potential energy surface of the reaction of the enzyme–inhibitor complex in which the Asp 25’ residue and N1 centre (pyrrolidine) are protonated while N2 centre (aziridine) remains deprotonated. The first reaction coordinate is the distance between the attacking carboxylate oxygen atom OD2 and the attacked C2 centre of the aziridine ring (distance A in Figure 9). The second reaction coChemPhysChem 2014, 15, 3226 – 3235

3232

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

The potential energy surface contains three distinct minima. The energetically highest minimum [E-Asp 25 ···I + H-Asp 25’] represents the reactant state of the reaction in which the ringopening reaction has not yet taken place (large value for the distance A) and the nitrogen centre of the aziridine ring is still deprotonated (large value for distance B). For the proton transfer from Aps 25’ to the aziridine nitrogen centre a very small barrier of only about 25 kJ mol 1 has to be surmounted. Due to the low barrier and because the resulting complex [EAsp 25 ···IH + + Asp 25’ ] is 18 kJ mol 1 more stable, one can assume that at the beginning of the reaction the proton is transferred from the Asp 25’ moiety to the aziridine ring (see Figure S8). For the subsequent ring opening reaction of 1020 (decreasing value for distance A) we Figure 9. QM/MM system of 1020 in its doubly protonated form. A and B indicate internal coordinates used to compute a reaction barrier of model the reaction coordinate. about 65 kJ mol 1. The geometrical structure of the transition state can be taken from Figure S9. It is higher than the barrier predicted for the direct ring-opening reaction of AZNP (20 kJ mol 1) but it is similar to the value predicted for the watermediated reaction of AZNP (50 kJ mol 1). The calculated reaction energy of about 30 kJ mol 1 indicates that this ring opening reaction of 1020 is considerably less exothermic than the corresponding reaction for AZNP ( 150 kJ mol 1). The geometrical structure of the covalent 1020 HIV-1 PR complex, that is, the final product of the inhibition reaction can be taken from Figure S10. Due to its pKa value of 7.9, also the nitrogen centre of the aziridine ring may be protonated under physiological conditions. If 1020 enters the active site doubly protonated, the preceding proton transfer from Asp 25’ cannot take place anymore. However, this is not necessary since the inhibition reaction diFigure 10. Contour plot of the PES for the bond formation between 1020 and the target enzyme. The definition of the distances A and B is given in Figure 9. Energies are given rectly starts from the position of the doubly protonin kJ mol 1. ated inhibitor [E-Asp 25 ···IH + ]. The actual ring-opening reaction corresponds to the second reaction step in Figure 10, for which a reaction energy of 32 kJ mol 1 and ordinate is B, which gives the distance between the nitrogen centre of the aziridine ring and the proton of the next water a barrier of 64 kJ mol 1 are predicted. MD simulations indicate molecule (Figure 9). To avoid artificial formations of energetithat a doubly protonated inhibitor fits into the active site. cally unfavorable oxonium or hydroxyl ions, the bond distanOur calculations predict that the second step of the inhibices B’ were adjusted at each point of the surface to allow the tion mechanism of 1020 has a reaction energy of only relaxation to the next lower local minimum. 30 kJ mol 1. Hence, the reaction is computed to be consider 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

ChemPhysChem 2014, 15, 3226 – 3235

3233

CHEMPHYSCHEM ARTICLES ably less exothermic than the corresponding reactions of known irreversible inhibitors (e.g. E64c/cathepsin B or EPNP/ HIV-1 PR). Their reaction energies were computed to be far below 120 kJ mol 1. This indicates that for 1020 the back reaction can take place, that is, 1020 is predicted to represent a covalent reversible inhibitor. Such a covalent reversible inhibitor should possess a similar efficacy as the corresponding covalent irreversible agent, but immunogenic side effects should be avoided.

3. Summary and Conclusions We present a protocol that allows a rational design of covalent inhibitors. The key structure in our approach is the non-covalent enzyme inhibitor complex, which is formed in the first step of the inhibition mechanism. It is important because its lifetime and geometry mainly influence the overall inhibition potency of the inhibitor. For example, the reaction leading to the covalent inhibitor–enzyme bond (second step of the inhibition mechanism) can only happen if the life-time of the noncovalent enzyme inhibitor complex is sufficiently long. The geometry of the non-covalent inhibitor complex is also important because long distances between the reacting centres of the enzyme and the inhibitor also diminish the reaction probability. Finally, the ability of the inhibitor to form such complexes with various enzymes is important for the selectivity of the inhibitor. Hence, in the first step of our protocol we compute the geometrical structure of the non-covalent enzyme inhibitor complex by means of QM/MM calculations. The calculations are necessary because the short lifetimes of these complexes prevent experimental characterizations. These short lifetimes result from the reactivity of the warheads. To circumvent this problem, our calculations start with the crystal structure of an appropriate covalent enzyme–inhibitor complex. The necessary structural information about the corresponding non-covalent enzyme–inhibitor complex is obtained by computing the second step of the inhibition mechanism backwards from the covalent to the non-covalent enzyme–inhibitor complex. After computing the structure of the non-covalent enzyme inhibitor complex of a given inhibitor, this structure is used as the starting point of docking experiments. These have the goal to obtain inhibitors with an improved affinity. Additionally, the warhead can also be changed either to obtain a better affinity or to increase the efficiency of subsequent bond formation step. Subsequently, using MD and QM/MM it is tested if the reaction is still possible for the new inhibitor. Using the information of the first round, further improvements can be achieved iteratively. We explored the possibilities of this protocol using EPNP as inhibitor for the HIV-1 PR. In our investigations we changed both, the recognition unit and the warhead of EPNP and could predict a new compound class which should show a much higher affinity within the non-covalent enzyme inhibitor complex. We used 1020 as an example for a potentially improved inhibitor to investigate the reactivity against the target protein. QM/MM calculations indeed indicate a sufficiently large reactivity. According to the calculated reaction energy 1020 is a cova 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org lent reversible inhibitor. Our protocol is an addition to available covalent docking procedures which allow retrieving potent covalent inhibitors of a given enzyme out of large compound libraries. To prove, refine or disprove our computational studies, experimental studies are under way.

Acknowledgements The authors thank the DFG for financial support via the SFB 630. P.K.S. acknowledges the Alexander von Humboldt Stiftung for valuable support. Keywords: covalent inhibition · drug discovery · hybrid approach · QM/MM · reaction mechanisms [1] I. D. Kuntz, K. Chen, K. A. Sharp, P. A. Kollman, Proc. Natl. Acad. Sci. USA 1999, 96, 9997 – 10002. [2] G. R. Desiraju, Acc. Chem. Res. 2002, 35, 565 – 573. [3] B. B. Brodie, W. D. Reid, A. K. Cho, G. Sipes, G. Krishna, J. R. Gillette, Proc. Natl. Acad. Sci. USA 1971, 68, 160 – 164. [4] J. Mitchell, D. Jollow, W. Potter, J. Gillette, B. Brodie, J. Pharmacol. Exp. Ther. 1973, 187, 211 – 217. [5] D. Jollow, J. Mitchell, W. Potter, D. Davis, J. Gillette, B. Brodie, J. Pharmacol. Exp. Ther. 1973, 187, 195 – 202. [6] W. Potter, D. Davis, J. Mitchell, D. Jollow, J. Gillette, B. Brodie, J. Pharmacol. Exp. Ther. 1973, 187, 203 – 210. [7] J. Mitchell, D. Jollow, W. Potter, D. Davis, J. Gillette, B. Brodie, J. Pharmacol. Exp. Ther. 1973, 187, 185 – 194. [8] A. J. T. Smith, X. Zhang, A. G. Leach, K. N. Houk, J. Med. Chem. 2009, 52, 225 – 233. [9] G. Jones, P. Willett, R. C. Glen, J. Mol. Biol. 1995, 245, 43 – 53. [10] G. Jones, P. Willett, R. C. Glen, A. R. Leach, R. Taylor, J. Mol. Biol. 1997, 267, 727 – 748. [11] M. Rarey, B. Kramer, T. Lengauer, G. Klebe, J. Mol. Biol. 1996, 261, 470 – 489. [12] C. Mcmartin, R. S. Bohacek, J. Comput. Aid. Mol. Des. 1997, 11, 333 – 344. [13] J. Schrçder, A. Klinger, F. Oellien, R. J. Marhçfer, M. Duszenko, P. M. Selzer, J. Med. Chem. 2013, 56, 1478 – 1490. [14] X. Ouyang, S. Zhou, C. T. T. Su, Z. Ge, R. Li, C. K. Kwoh, J. Comput. Chem. 2013, 34, 326 – 336. [15] G. M. Morris, R. Huey, W. Lindstrom, M. F. Sanner, R. K. Belew, D. S. Goodsell, A. J. Olson, J. Comput. Chem. 2009, 30, 2785 – 2791. [16] M. Mladenovic, K. Ansorg, R. F. Fink, W. Thiel, T. Schirmeister, B. Engels, J. Phys. Chem. B 2008, 112, 11798 – 11808. [17] M. Mladenovic, K. Junold, R. F. Fink, W. Thiel, T. Schirmeister, B. Engels, J. Phys. Chem. B 2008, 112, 5458 – 5469. [18] J. Tang, J. Biol. Chem. 1971, 246, 4510 – 4517. [19] T. D. Meek, B. D. Dayton, B. W. Metcalf, G. B. Dreyer, J. E. Strickler, J. G. Gorniak, M. Rosenberg, M. L. Moore, V. W. Magaard, C. Debouck, Proc. Natl. Acad. Sci. USA 1989, 86, 1841 – 1845. [20] R. Salto, L. M. Bab, J. Li, J. R. Ros, Z. Yu, A. Burlingame, J. J. D. Voss, Z. Sui, P. O. d. Montellano, C. S. Craik, J. Biol. Chem. 1994, 269, 10691 – 10698. [21] T. Schmidt, A. Paasche, C. Grebner, K. Ansorg, J. Becker, W. Lee, B. Engels, Top. Curr. Chem. 2012, 1 – 77. [22] H. M. Senn, W. Thiel, Angew. Chem. Int. Ed. 2009, 48, 1198 – 1229; Angew. Chem. 2009, 121, 1220 – 1254. [23] M. Mladenovic, T. Schirmeister, S. Thiel, W. Thiel, B. Engels, ChemMedChem 2007, 2, 120 – 128. [24] M. Mladenovic, R. F. Fink, W. Thiel, T. Schirmeister, B. Engels, J. Am. Chem. Soc. 2008, 130, 8696 – 8705. [25] U. Ryde, L. Olsen, K. Nilsson, J. Comput. Chem. 2002, 23, 1058 – 1070. [26] U. Ryde, K. Nilsson, J. Am. Chem. Soc. 2003, 125, 14232 – 14233. [27] U. Ryde, K. Nilsson, J. Mol. Struct. THEOCHEM 2003, 632, 259 – 275. [28] V. Buback, M. Mladenovic, B. Engels, T. Schirmeister, J. Phys. Chem. B 2009, 113, 5282 – 5289.

ChemPhysChem 2014, 15, 3226 – 3235

3234

CHEMPHYSCHEM ARTICLES [29] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kal, K. Schulten, J. Comput. Chem. 2005, 26, 1781 – 1802. [30] B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, M. Karplus, J. Comput. Chem. 1983, 4, 187 – 217. [31] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, M. L. Klein, J. Chem. Phys. 1983, 79, 926 – 935. [32] V. Zoete, M. A. Cuendet, A. Grosdidier, O. Michielin, J. Comput. Chem. 2011, 32, 2359 – 2368. [33] K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov, A. D. MacKerell, J. Comput. Chem. 2010, 31, 671 – 690. [34] H. K. Hall, J. Am. Chem. Soc. 1957, 79, 5441 – 5444. [35] R. K. Bansal, Heterocyclic Chemistry, Anshan Ltd., TunbridgeWells, 2008. [36] ChemShell, a computational chemistry shell, see http://www.chemshell.org. [37] P. Sherwood, A. H. de Vries, M. F. Guest, G. Schreckenbach, C. A. Catlow, S. A. French, A. A. Sokol, S. T. Bromley, W. Thiel, A. J. Turner, S. Billeter, F. Terstegen, S. Thiel, J. Kendrick, S. C. Rogers, J. Casci, M. Watson, F. King, E. Karlsen, M. Sjovoll, A. Fahmi, A. Schaefer, C. Lennartz, J. Mol. Struct. THEOCHEM 2003, 632, 1 – 28. [38] J. Kstner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander, P. Sherwood, J. Phys. Chem. A 2009, 113, 11856 – 11865. [39] TURBOMOLE v6.3.1, a development of university of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989 – 2013, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. [40] A. D. Becke, Phys. Rev. A 1988, 38, 3098 – 3100. [41] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785 – 789. [42] A. Schfer, C. Huber, R. Ahlrichs, J. Chem. Phys. 1994, 100, 5829 – 5835. [43] M. L. Verdonk, J. C. Cole, M. J. Hartshorn, C. W. Murray, R. D. Taylor, Proteins 2003, 52, 609 – 623. [44] C. A. Sotriffer, P. Sanschagrin, H. Matter, G. Klebe, Proteins 2008, 73, 395 – 419. [45] R. B. Rose, C. S. Craik, N. L. Douglas, R. M. Stroud, Biochemistry 1996, 35, 12933 – 12944. [46] R. B. Rose, J. R. Rose, R. Salto, C. S. Craik, R. M. Stroud, Biochemistry 1993, 32, 12498 – 12507. [47] E. Rutenber, E. B. Fauman, R. J. Keenan, S. Fong, P. S. Furth, P. R. O. d. Montellano, E. Meng, I. D. Kuntz, D. L. DeCamp, R. Salto, J. Biol. Chem. 1993, 268, 15343 – 15346. [48] E. E. Rutenber, F. McPhee, A. P. Kaplan, S. L. Gallion, J. C. Hogan Jr., C. S. Craik, R. M. Stroud, Bioorg. Med. Chem. 1996, 4, 1545 – 1558. [49] H. O. Andersson, K. Fridborg, S. Loewgren, M. Alterman, A. Mhlman, M. Bjoersne, N. Garg, I. Kvarnstrm, W. Schaal, B. Classon, A. Karln, U. H. Danielsson, G. Ahlsn, U. Nillroth, L. Vrang, B. Oeberg, B. Samuelsson, A. Hallberg, T. Unge, Eur. J. Biochem. 2003, 270, 1746 – 1758. [50] T. Sklov, J. Hasˇek, J. Dohnlek, H. Petrokov, E. Buchtelov, J. Dusˇkov, M. Soucˇek, P. Majer, T. Uhlkov, J. Konvalinka, J. Med. Chem. 2003, 46, 1636 – 1644. [51] J. K. Ekegren, N. Ginman, A. Johansson, H. Wallberg, M. Larhed, B. Samuelsson, T. Unge, A. Hallberg, J. Med. Chem. 2006, 49, 1828 – 1832. [52] H. E. Klei, K. Kish, P.-F. M. Lin, Q. Guo, J. Friborg, R. E. Rose, Y. Zhang, V. Goldfarb, D. R. Langley, M. Wittekind, S. Sheriff, J. Virol. 2007, 81, 9525 – 9535.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org [53] Y. Tie, P. I. Boross, Y.-F. Wang, L. Gaddis, A. K. Hussain, S. Leshchenko, A. K. Ghosh, J. M. Louis, R. W. Harrison, I. T. Weber, J. Mol. Biol. 2004, 338, 341 – 352. [54] S. Chellappan, G. S. Kiran, K. Reddy, A. Ali, M. N. L. Nalam, S. G. Anjum, H. Cao, V. Kairys, M. X. Fernandes, M. D. Altman, B. Tidor, T. M. Rana, C. A. Schiffer, M. K. Gilson, Chem. Biol. Drug Des. 2007, 69, 298 – 313. [55] S. Thaisrivongs, K. D. Watenpaugh, W. J. Howe, P. K. Tomich, L. A. Dolak, K.-T. Chong, C.-S. C. Tomich, A. G. Tomasselli, S. R. Turner, J. Med. Chem. 1995, 38, 3624 – 3637. [56] X. Wu, P. Oehrngren, J. K. Ekegren, J. Unge, T. Unge, H. Wallberg, B. Samuelsson, A. Hallberg, M. Larhed, J. Med. Chem. 2008, 51, 1053 – 1057. [57] A. Mahalingam, L. Axelsson, J. K. Ekegren, J. Wannberg, J. Kihlstrom, T. Unge, H. Wallberg, B. Samuelsson, M. Larhed, A. Hallberg, J. Med. Chem. 2010, 53, 607 – 615. [58] J. D. A. Tyndall, L. K. Pattenden, R. C. Reid, S.-H. Hu, D. Alewood, P. F. Alewood, T. Walsh, D. P. Fairlie, J. L. Martin, Biochemistry 2008, 47, 3736 – 3744. [59] N. M. King, M. Prabu-Jeyabalan, R. M. Bandaranayake, M. N. L. Nalam, E. A. Nalivaika, A. zen, T. Halilolu, N. K. Ylmaz, C. A. Schiffer, ACS Chem. Biol. 2012, 7, 1536 – 1546. [60] M. N. L. Nalam, A. Ali, M. D. Altman, G. S. K. K. Reddy, S. Chellappan, V. Kairys, A. zen, H. Cao, M. K. Gilson, B. Tidor, T. M. Rana, C. A. Schiffer, J. Virol. 2010, 84, 5368 – 5378. [61] A. K. Ghosh, S. Kulkarni, D. D. Anderson, L. Hong, A. Baldridge, Y.-F. Wang, A. A. Chumanevich, A. Y. Kovalevsky, Y. Tojo, M. Amano, Y. Koh, J. Tang, I. T. Weber, H. Mitsuya, J. Med. Chem. 2009, 52, 7689 – 7705. [62] R. Ishima, Q. Gong, Y. Tie, I. T. Weber, J. M. Louis, Proteins Struct. Funct. Bioinf. 2010, 78, 1015 – 1025. [63] F. M. Olajuyigbe, N. Demitri, J. O. Ajele, E. Maurizio, L. Randaccio, S. Geremia, Med. Chem. Lett. 2010, 1, 254 – 257. [64] Y. Kawasaki, E. E. Chufan, V. Lafont, K. Hidaka, Y. Kiso, L. M. Amzel, E. Freire, Chem. Biol. Drug Des. 2010, 75, 143 – 151. [65] C.-H. Shen, Y.-F. Wang, A. Y. Kovalevsky, R. W. Harrison, I. T. Weber, FEBS J. 2010, 277, 3699 – 3714. [66] M. A. Navia, P. M. D. Fitzgerald, B. M. McKeever, C.-T. Leu, J. C. Heimbach, W. K. Herber, I. S. Sigal, P. L. Darke, J. P. Springer, Nature 1989, 337, 615 – 620. [67] A. Wlodawer, M. Miller, M. Jaskolski, B. K. Sathyanarayana, E. Baldwin, I. T. Weber, L. M. Selk, L. Clawson, J. Schneider, S. B. Kent, Science 1989, 245, 616 – 621. [68] V. Hornak, C. Simmerling, Drug Discovery Today 2007, 12, 132 – 138. [69] S. R. Turner, J. W. Strohbach, R. A. Tommasi, P. A. Aristoff, P. D. Johnson, H. I. Skulnick, L. A. Dolak, E. P. Seest, P. K. Tomich, M. J. Bohanon, M.-M. Horng, J. C. Lynn, K.-T. Chong, R. R. Hinshaw, K. D. Watenpaugh, M. N. Janakiraman, S. Thaisrivongs, J. Med. Chem. 1998, 41, 3467 – 3476. [70] M. J. Turner, J. J. McKinnon, D. Jayatilaka, M. A. Spackman, CrystEngComm 2011, 13, 1804 – 1813.

Received: July 25, 2014 Revised: August 29, 2014 Published online on September 24, 2014

ChemPhysChem 2014, 15, 3226 – 3235

3235