Andrew G.
Leach
*a,
Lise-Lotte
Olsson
b and
Daniel J.
Warner
c
aAstraZeneca, 30S373 Mereside, Alderley Park, Macclesfield, SK10 4TG, UK. E-mail: andrew.leach@astrazeneca.com; Tel: +44 (0)1625 231853
bAstraZeneca, Pepparedsleden 1, 431 83 Mölndal, Sweden. E-mail: lise-lotte.olsson@astrazeneca.com; Tel: +46 31 7064931
cAstraZeneca, 3328, 7171 Frederick Banting, Ville St. Laurent, Montreal, Canada
First published on 14th September 2012
Two series of compounds that inhibit dimerisation of iNOS have been reported in the literature and studied crystallographically. We have applied a range of docking techniques to these compounds and shown that the geometry of the complexes can be reproduced well but that even within a series of related compounds the docking scores do not rank compounds well. Quantum mechanical studies using a model system for the protein and a range of density functionals are able to reproduce the geometry of the complexes and the corresponding affinities. The combination of docking to generate geometries and quantum mechanical calculations to predict complexation energies is a powerful one for rationalising observed changes and designing improved compounds.
A number of compounds have been described in the literature as inhibitors of iNOS but these have generally been mimics of arginine, the substrate for the enzyme, which is present at high physiological concentrations.2 Such molecules have not shown a suitable balance of efficacy and tolerability for clinical progression. This may be because they are unable to compete with the high concentration of natural substrate and may bind at other sites that recognise arginine.
This background has spurred an interest in compounds that inhibit iNOS in ways which are not competitive with arginine. Two series of compounds have been disclosed by workers from Berlex and Kalypsys which were described as inhibiting the dimerisation of iNOS.3–5 Dimerisation is a prerequisite for enzyme activity and the process requires the cofactor tetrahydrobiopterin which also delivers electrons to the heme iron that is present in the enzyme active site.6 This initiates the reaction that leads first to binding of oxygen and then to oxidation of arginine to give NO and citrulline.1
The compounds described by Berlex derived from a screening library designed to include a group that binds to the heme iron. Two representatives of this series had been studied with protein–ligand crystal structures in which the protein had been engineered to not dimerise.3,7 The structures of these compounds (1 and 2) in complex with murine iNOS are illustrated in Fig. 1. This shows that a principal feature of these compounds' binding mode is a bond between the heme iron of iNOS and an imidazole nitrogen. An intriguing further feature is the positioning of the dioxolane group in both compounds such that one oxygen is able to form a weak interaction with a backbone NH and one of the polarized CH bonds a weak interaction with a backbone carbonyl. These hydrogen bonding features are shared with the binding mode of compound 3 described below.
Fig. 1 Protein–ligand crystal structure of compounds 1 and 2 in complex with monomeric iNOS.8 |
The compounds described by Kalypsys derived from high throughput screening.4 Building on the work that yielded the protein–ligand crystal structures of 1 and 2, we were able to obtain the protein–ligand crystal structure of an exemplar of the series of compounds described by Kalypsys, compound 3 (Fig. 2).5 By contrast to 1 and 2, compound 3 reported by Payne et al. is found not to form an interaction with the heme iron.5 It relies upon hydrogen bonding between the carbonyl of 3 and the backbone NH of Met368 and the amidic NH of 3 and the backbone carbonyl of Trp366. These interactions are analogous to those involving the weak donor and acceptor of the dioxolane group in 1 and 2. A further key aspect of the interaction between 3 and iNOS is that the aromatic core of 3 stacks onto the extended aromatic system of the porphyrin in the heme group. The contribution of such interactions to binding energies is challenging to quantitate. The imidazopyrazine fragment of 3 protrudes into a water filled part of the pocket and does not make any direct interactions with the protein.
Fig. 2 Protein–ligand crystal structure of compound 3 in complex with monomeric iNOS.8 |
The assortment of interaction types involved in the interaction between iNOS monomer and the inhibitors described by both Berlex and Kalypsys make this an extremely complex system for which to make predictions about binding pose and binding energy. It serves as a good test to identify protocols that are able to treat the balance between metal–ligand interactions, hydrogen bonds and π-stacking correctly. There are two stages to achieving this aim, the first is to find a method that predicts the binding pose correctly, the second is to find a method that predicts the binding energy correctly. The first has been approached by studying a small range of commonly available docking protocols. When considering approaches to the second problem, it might be expected that quantum mechanical calculations can achieve an adequate balance between the different interaction types. However, there are many challenges to doing this. An appropriate model system to probe the SAR must be identified and an appropriate level of theory must be found. Hartree–Fock calculations are likely to perform poorly on this balance of interaction types and many of the density functional theory approaches are known to treat dispersion based interactions poorly.9 Recently, the Minnesota set of functionals have offered improvements in this area10,11 and the ωB97XD functional includes an empirical correction to deal with dispersion.12 Beyond these considerations it has also been of interest to investigate whether the computational methods used are able to describe the SAR for either IC50 or lipophilicity ligand efficiency. This latter property takes into account the enhanced binding to proteins that can come about by simply increasing the lipophilicty of a molecule and hence its tendency to partition into more hydrophobic environments and away from water.
Fig. 3 A series of compounds described by Kalypsys which all share a common sidechain.4 Quantum mechanical modelling of these compounds abbreviated this sidechain to methyl and this is designated with (′). |
The geometry obtained for compounds 1 to 3 by each method were examined and compared to their crystal structures. For compounds 1 and 2, three distances were measured: the Fe–N distance, the distance between the dioxolane oxygen and the backbone nitrogen of Met368 and the distance between the dioxolane carbon and the backbone carbonyl oxygen of Trp366 (the last two interaction distances are indicated in Fig. 1). For compound 3 the distances between the ligand carbonyl and the backbone nitrogen of Met368 and between the ligand nitrogen (of its NH functional group) and the backbone carbonyl oxygen of Trp366 (the two distances indicated in Fig. 2) were measured in the docked poses. The RMSDs for two fragments of each ligand were also measured after the crystal structures of 1 and 2 and the docked poses for all structures were aligned to that for 3 using the inner (16-membered) ring of the porphyrin. For 1 and 2, the RMSDs for the imidazole-pyrimidine and the benzodioxole estimate how far away these fragments are from where they are in the crystal structures while for 3 the RMSDs for the 10 ring atoms of the quinolinone and for the 13 atoms of the cyclobutylsubstituted imidazopyrazine were computed. All of these parameters are shown in Table S1 in the ESI.†
The geometric parameters in Table S1† show that the poses that GOLD ranks best for either 1 or 2 are incorrect; they do not include any interaction with the heme iron. All of the GOLD protocols are however able to locate a pose for 1 that involves an interaction between the imidazole N and the heme iron when docking into either of P1 or P2, but do not rank these poses best. GOLD is only able to find a pose for compound 2 that includes an interaction with the metal when docking into its own protein P2 and using the GOLD3 method. By contrast, the Glide protocol XP1 is able to correctly position each ligand in its own crystal structure. None of the methods are able to reliably position any of the ligands in any of the proteins other than that from which it was extracted. Examination of the protein structures does not reveal any drastic differences between them other than Glu371 in the complex of 3 which occludes the benzodioxole binding region of 1 and 2 (Fig. 4). The heme also varies slightly both in the orientation of its sidechains and the degree of twisting; the iron atom is depressed away from the ligand in the complex of 3 presumably because it is not forming a direct interaction with the ligand. It remains an impressive achievement that one of the methods is able to correctly place this set of ligands which utilise such a complex array of interaction types when binding to iNOS in its monomeric form.
Fig. 4 Protein–ligand complex structures of compounds 1 (orange), 2 (brown) and 3 (blue) overlaid using the inner ring of the heme porphyrin. |
As the docking protocols are able to dock compound 3 well, albeit with some ambiguity concerning the orientation of the imidazopyrazine, it seemed reasonable that they should also be able to position compounds 4–16 which are expected to bind with the quinolinone in the same position as in 3. This is borne out in practice and the poses for most compounds are those expected although it was notable that only one of the Glide protocols was able to dock the methylated quinolinone 9. The scores arising from each pose are presented in Table S2 in the ESI† and summarised in Table 1 and show that although the binding poses are similar to that observed for 3, the resulting scores do not correlate well with either pIC50 or LLE. The best correlation for pIC50 is with Glide docking into P2 using XP1 which gives an R2 value of 0.43 and RMSE of 1.40. For LLE, docking into P2 using XP1 gives R2 of 0.6 and RMSE of about 1.8. When the experimental values that fall outside the range of the assay are excluded from the correlation (compounds 5, 8, 9 and 10), the best correlation for pIC50 is for Glide docking into P3 with SP1 which method gives R2 of 0.68 and RMSE of 0.77. For LLE it is XP1 docking into P2 which gives R2 of 0.61 and RMSE of 1.44, although SP1 using P3 gives similar performance. Given the observations above concerning docking poses, it is of note that this dataset is dominated by compounds that would be expected to dock most similarly to compound 3; it is docking into that structure which is most likely to give the right correlation for the right reasons.
Protein | Method | pIC50 | LLE | ||
---|---|---|---|---|---|
R 2 | RMSE | R 2 | RMSE | ||
P1 | SP1 | 0.04 | 1.81 | 0.08 | 2.67 |
XP1 | 0.01 | 1.84 | 0.02 | 2.76 | |
SP2 | 0.22 | 1.63 | 0.10 | 2.64 | |
XP2 | 0.00 | 1.85 | 0.06 | 2.71 | |
GOLD1 | 0.26 | 1.64 | 0.26 | 2.42 | |
GOLD2 | 0.25 | 1.65 | 0.45 | 2.08 | |
GOLD3 | 0.20 | 1.71 | 0.17 | 2.56 | |
P2 | SP1 | 0.08 | 1.77 | 0.04 | 2.73 |
XP1 | 0.43 | 1.40 | 0.60 | 1.77 | |
SP2 | 0.30 | 1.55 | 0.35 | 2.24 | |
XP2 | 0.20 | 1.70 | 0.13 | 2.63 | |
GOLD1 | 0.06 | 1.85 | 0.05 | 2.75 | |
GOLD2 | 0.24 | 1.66 | 0.41 | 2.16 | |
GOLD3 | 0.18 | 1.73 | 0.15 | 2.60 | |
P3 | SP1 | 0.29 | 1.56 | 0.32 | 2.31 |
XP1 | 0.06 | 1.79 | 0.15 | 2.57 | |
SP2 | 0.01 | 1.84 | 0.00 | 2.79 | |
XP2 | 0.01 | 1.84 | 0.02 | 2.77 | |
GOLD1 | 0.13 | 1.77 | 0.01 | 2.80 | |
GOLD2 | 0.13 | 1.78 | 0.25 | 2.44 | |
GOLD3 | 0.06 | 1.85 | 0.00 | 2.82 |
To perform quantum mechanical calculations on a system such as the iNOS monomer, it is necessary to extract a model for the binding site that includes the atoms that are most important for interacting with the ligands involved. In particular, this has focused on the set of ligands 4 to 16. The binding pose of these compounds is expected to be close to that of 3 which suggested that a good description of the binding site could be provided by: (i) the backbone of residues 364–368 (ii) the Cys that interacts with the iron reduced to SMe (iii) iron porphyrin with the various sidechains on the heme porphyrin ligand removed. The required atoms were obtained by editing the protein–ligand crystal structure of 3 and are shown in Fig. 5. In order to mimic the scaffolding effect of the rest of the protein, the carbon and nitrogen atoms of the back bone (coloured orange and purple respectively) as well as the four bridging carbons of the porphyrin which are chosen to be remote from the iron centre both in terms of distance and electronically are all fixed in space using the fixed Cartesian redundant coordinate feature in Gaussian09. To further simplify these studies, the common sidechain to compounds 4–16 was abbreviated to a methyl group to give compounds 4′ to 16′. For these compounds because the rest of the molecule remains the same, the methyl substituted compounds should have parallel SAR to the fully elaborated compounds. Compounds 1 and 2 that interact with the heme iron directly and the backbone of the protein through parts of the molecule that are remote from one another were studied as the full molecule. Initial positioning of these two compounds was performed by overlaying the heme inner ring of their protein–ligand complexes with that of the QM model system. It is expected that because the model receptor excludes Glu371 and allows the porphyrin ligand to change shape that it should allow 1 and 2 to bind appropriately.
Fig. 5 The quantum mechanical model system extracted from the complex of 3 with monomeric iNOS. |
One matter that crystal structures are not able to probe directly is the spin state of iron. While it is assumed in this case that the iron is in its +3 oxidation state, it remains possible for it to adopt a high spin sextet or low spin doublet electronic configuration. It may also be able to adopt an intermediate spin quartet but this possibility is more difficult to treat computationally and has not been considered further. In order to address this issue, the calculations have been performed for both spin states of the iron.
A range of quantum mechanical methods were used in these studies. A set of calculations using HF were performed although these had to be combined with the quadratic convergence protocol which slows the calculations down. The B3LYP functional has long been favoured as showing a good balance of speed and accuracy although it is known to treat dispersion based interactions poorly.9,22,23 The M06, M06-2X, M06L and M06HF functionals have all been described by the Truhlar group.10,11 M06 is described as having general applicability, M06-2X as being good at treating aromatic–aromatic interactions but poor at transition metal chemistry, M06L as more accurate for transition metals and M06HF for a number of weakly bonded states. Because it includes an HF component M06HF was also used with the quadratic convergence protocol. The iNOS system represents a good test for this set of functionals. Finally, the ωB97XD includes an empirical correction for dispersion.12 All of these functionals were combined with the 6-31G* basis set for all atoms apart from iron which used the lanl2dz basis set and pseudopotential.24–26 All calculations were performed in Gaussian09 and ligand geometries were also optimized and characterized by frequency calculations whereas calculations involving the Cartesian constraints did not include frequency calculations.27 The geometry of all complexes and ligands obtained by these calculations are provided as ESI.†
The geometry of the complexes of compounds 1 to 3 (the latter modelled by 13′) obtained by each method was examined and the same distances as studied for the docking methods were measured. In addition the Fe–S distance (which was not variable in the docking studies) which depends on the spin state of the iron was measured. These distances are reported in Table S3 in the ESI† and show that for all methods and all three compounds the distances found in the QM system for the doublet state of iron agree better with those observed in the crystal structures and suggest that in the crystallographic construct the iron is in the low spin state regardless of whether the ligand forms a direct interaction with the metal or not. Of the various methods studied, the M06 and M06L methods fare best overall in terms of predicting the bond lengths studied whereas M06-2X fares best for the set of distances excluding Fe–S. M06-2X is the only method that predicts that the Fe–S distance should be shorter in the sextet than the doublet and also has a smaller difference between the spin states for the Fe–N distance than other methods, this supports the suggestion that this method is less good for transition metal systems.
Finally the energy of complexation for each ligand binding to the QM model system in the two spin states has been computed, using separate optimizations of an empty receptor with the Cartesian constraints still applied. Full optimization of the isolated ligands complete the set of components needed to compute the energy change, Ecomplex, that corresponds to eqn (1).
(1) |
The geometry of the complex for compounds 4′ to 16′ was generally obtained by editing that of the model for compound 13′ which in turn came from the crystal structure of compound 3. The ligand geometries were extracted from the complex but as with the docking, corina was used to generate a conformation for compounds 1 and 2. In both cases following optimization by all methods the geometry obtained from the complex was found to be lower in energy than that generated by corina. Compound 6′ has two possible chair conformations and both were generated for the isolated ligand and the complex and the lowest energy conformation of each was used in the resulting energy changes. The other model compounds do not have any conformational degrees of freedom and so were not studied in any more detail.
The complexation energies for the different ligands to the model receptor in doublet and sextet states are provided in Table S4 in the ESI† and performance statistics summarise this in Table 2.28 In general the various methods perform well, with the best method for pIC50 being M06HF with the receptor in the doublet state. This method gives an R2 of 0.67 and RMSE of 0.96. By comparison, for the same set of compounds, the best docking method is GOLD1 docking into P3 which gives R2 of 0.67 and RMSE of 0.98. When LLE is considered instead, the quantum mechanical methods perform significantly better and the best is once again M06HF binding to the doublet state which gives R2 of 0.94 and RMSE of 0.46; the accuracy is similar to that expected in most biological assays. Many of the QM methods give excellent R2 values over 0.8 and RMSE below 0.9. By contrast for the same set of LLE values, the best docking method is GOLD3 docking into P3 which gives R2 of only 0.47 and RMSE of 1.56. The QM approach is sufficiently accurate to guide the design of compounds in lead optimisation where small structural changes leading to small changes in binding energy are attempted and where docking methods are often found wanting. It is worth considering why the QM methods perform so much better when modelling LLE than pIC50. The model system for the QM calculations is actually each of the species shown in eqn (1) in the gas phase. No account is taken of the solvation energy of the different ligands. The clogP term takes this into account to a certain degree, to the extent that 1-octanol provides a good model for the gas phase.
Methods | Spin state | pIC50 | LLE | ||
---|---|---|---|---|---|
R 2 | RMSE | R 2 | RMSE | ||
B3LYP | Doublet | 0.19 | 1.53 | 0.52 | 1.48 |
HF | Doublet | 0.28 | 1.46 | 0.41 | 1.65 |
M06 | Doublet | 0.53 | 1.16 | 0.79 | 0.98 |
M062X | Doublet | 0.64 | 1.00 | 0.84 | 0.85 |
M06HF | Doublet | 0.67 | 0.96 | 0.94 | 0.46 |
M06L | Doublet | 0.66 | 0.99 | 0.81 | 0.94 |
ωB97XD | Doublet | 0.61 | 1.06 | 0.78 | 0.99 |
B3LYP | Sextet | 0.22 | 1.51 | 0.52 | 1.48 |
M06 | Sextet | 0.50 | 1.20 | 0.76 | 1.04 |
M062X | Sextet | 0.63 | 1.03 | 0.84 | 0.86 |
M06L | Sextet | 0.61 | 1.06 | 0.79 | 0.99 |
ωB97XD | Sextet | 0.27 | 1.45 | 0.37 | 1.69 |
Although it is gratifying to find computational methods that are able to correctly rank compounds, it is arguably of more use to obtain structural insight from those calculations to guide future design. For instance, there is an increase in LLE upon addition of an extra fluorine atom going from 4 to 13 and also on going from 12 to 16. This is mirrored in an increased Ecomplex computed by M06HF. Examination of the geometries obtained by M06HF show that the more fluorinated compound is moved slightly away from the protein backbone (the carbon bearing H or F moves by 0.15 Å) suggesting weaker interaction at that point while the donor–acceptor pair move less which suggests that the increased Ecomplex must be due to an improvement in the π–π interaction caused by addition of fluorine.
It is of interest that the binding energies of compounds 1 and 2 are predicted to be smaller than those of model compounds 4′ to 16′ despite the fact that 1 and 2 are measured to be tight binding compounds. It is possible that this discrepancy arises because of errors introduced by reducing to model compounds, or an inadequacy of the description of the gap between doublet and sextet states or due to a kinetic barrier to spin flipping that prevents some of the compounds from achieving their lowest energy state (Scheme 1). Although the sextet state is generally higher in energy than the doublet, it is possible that only the sextet state is kinetically accessible. For compounds 4′ to 16′, the gap between the doublet and the sextet is computed to change by no more than about 1 kcal mol−1 upon complexation. However, for compounds 1 and 2, the gap is made larger by around 8–12 kcal mol−1 according to B3LYP, M06 and M06L. This would provide a significantly higher driving force for the spin flip for these two compounds than for the others. If the compounds that do not form direct bonds to iron (4 to 16) are not able to access the doublet state of the complex whereas compounds, like 1 and 2 can then this would provide them a significant extra binding energy. It may also instill these compounds with slow off rate kinetics. This interesting possibility could be the subject of experimental studies.
Scheme 1 Potential mechanism of inhibition for compounds that form a direct interaction with iron not available to other ligand types. |
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c2md20159c |
This journal is © The Royal Society of Chemistry 2013 |