Beyond Traditional Structure-Based Drug Design: The Role of Iron Complexation, Strain, and Water in the Binding of Inhibitors for Hypoxia-Inducible Factor Prolyl Hydroxylase 2.

A combination of structure-based drug design and medicinal chemistry efforts led us from benzimidazole-2-carboxamide with modestly active hypoxia-inducible factor prolyl hydroxylase 2 inhibition to certain benzimidazole-2-pyrazole carboxylic acids that were more potent as well as orally efficacious stimulators of erythropoietin secretion in our in vivo mouse model. To better understand the structure-activity relationship, it was necessary to account for (i) the complexation of the ligand with the active site Fe2+, (ii) the strain incurred by the ligand upon binding, and (iii) certain key water interactions identified by a crystal structure analysis. With this more complete computational model, we arrived at an overarching paradigm that accounted for the potency differences between benzimidazole-2-carboxamide and benzimidazole-2-pyrazole carboxylic acid enzyme inhibitors. Moreover, the computational paradigm allowed us to anticipate that the bioisostere replacement strategy (amide → pyrazole), which had shown success in the benzimidazole series, was not generally applicable to other series. This illustrates that to fully reconcile the important ligand-active site interactions for certain targets, one often needs to move beyond traditional structure-based drug design (such as crystallographic analysis, docking, etc.) and appeal to a higher level of computational theory.

Introduction

The hypoxia-inducible factor prolyl hydroxylase domain enzymes (HIF-PHD) are a family of highly conserved, iron-containing enzymes that play a critical role in adaptation and survival during oxygen-deficient conditions,[1−3] with the prolyl hydroxylase 2 (PHD2) isoform being of particular interest in drug discovery since it holds promise for the treatment of anemia, myocardial infarction, stroke, and metabolic disorders through inhibition of degradation of the HIF transcription factors.[4]

A combination of virtual screening, structure-based drug design, and medicinal chemistry efforts led to certain benzimidazole-2-carboxamide compounds, which showed modest enzyme inhibition[5] (e.g., 1, pIC50 = 4.9) and generally failed to significantly stimulate erythropoietin (EPO) release (a major transcriptional factor for HIF) from Hep3B cells at concentrations up to 100 μM. Indeed, overall optimization of the in vitro potency of the benzimidazole-2-carboxamide series proved challenging.

A bioisostere replacement[6] effort eventually led to the benzimidazole-2-pyrazole carboxylic acid series, which were significantly more potent compared to benzimidazole-2-carboxamide with several compounds having pIC50 ≥ 7.0. Moreover, many of these compounds increased both HIF-1α accumulation and EPO release in Hep3B cells in vitro,[7] and some produced an increase in EPO when orally dosed in vivo (mouse model).[8]

Despite these successes, we desired to have a clearer understanding of the underlying mechanisms responsible for the dramatic potency improvement upon bioisosteric replacement of the amide with pyrazole and whether this was indicative of a general strategy for all amide-based PHD2 enzyme inhibitors.

By systematically refining the computational model and level of theory, we were able to account for this improvement in enzyme potency and lay the foundation for a more general binding hypothesis.

Results and Discussion

Active Site Interactions of the Benzimidazole Series

To better understand the ligand–active site interactions, we obtained a crystal structure of 5,6-dichloro-benzimidazole-2-carboxamide (1) complexed with PHD2 (PDB ID: 3OUI). Although 1 is the most potent compound from this series, its pIC50 was still a modest value of 4.9. Figure depicts the key interactions of 1: (i) a salt bridge with R383 via the acid moiety, (ii) a bridging water interaction with Y303 via the NH of benzimidazole, and (iii) a bidentate interaction with Fe2+ via the nitrogen atom of benzimidazole and oxygen atom of the amide.

Figure 1

Crystal structure of 1 in the active site of PHD2 (PDB ID: 3OUI). In addition to the bidentate interaction with Fe2+, the other key interactions are shown as dashed lines.

Despite being the more potent series, members of the benzimidazole-2-pyrazole carboxylic acid series also revealed the same key interactions as 1 (e.g., 2c in Figure ) upon obtaining the crystal structures.

Figure 3

Crystal structure of 2c in the active site of PHD2 (PDB ID: 3OUH). In addition to the bidentate interaction with Fe2+, the other key interactions are shown as dashed lines. (Note that the same interactions are also present in 1 as seen in Figure .) The combination of docking and Fe2+ calculations correctly predicted this binding pose and orientation of the disubstitution prior to obtaining the crystal structure.

Thus, while the crystallography information provided valuable insight into the important interactions occurring in the benzimidazole series, neither it nor the subsequent docking analysis reconciled the two-log unit difference in pIC50, which resulted from replacing the amide motif with the pyrazole bioisostere. To better understand this, it was clear that a computational model that evolved past the usual structure-based approaches of docking and crystallographic analyses would be necessary. Any learning could then potentially be utilized to identify a new chemical series with increased potency.

While docking captures the general shape of the active site and key interactions at a force-field level of theory, we anticipated that several important binding characteristics would be missed. In particular, the Fe2+–ligand interaction is not captured at a high level of theory. Moreover, both docking and crystallography results of 1 and 2 revealed that a near-planar binding pose was necessary to facilitate the bidentate complexation of Fe2+, whereas our initial computational analysis of the unbound ligand (gas phase) revealed a nonplanar conformation; this implied that a strain penalty upon binding would be incurred, which warranted a more detailed computational analysis. Finally, we felt that the bridging water interaction and its contribution to the overall potency also necessitated further investigation.

Complexation and Strain

To more thoroughly interrogate our system, we supplemented our initial analyses with quantum calculations using Jaguar (under Maestro Version 10.7.015, release 2016-3)[9] with density functional theory (B3LYP) and the LACV3P* basis set. Further, symmetry was turned off, the grid density was set to “fine”, and the accuracy level was set to “accurate”. At this level of theory, we anticipated that we would gain deeper insight into the Fe2+–ligand complexation and binding strain. We initially explored a model system that would represent the key residues involved in complexation (H313, D315, and H374), Fe2+, the complexing water (known from crystallography), and the ligand. Here, the histidine residues were approximated as imidazoles and the aspartate residue as acetate. To further refine the model and since it was common to both benzimidazole analogs, we ignored the carboxylic acid of the ligands. This had the added benefit of simplifying our calculations by eliminating the need for diffuse functions (which can cause convergence issues).

Using this initial model system, we found that the total energy for both Fe2+ complexation and strain was lower (more favorable) for the pyrazole motif, and therefore, it was consistent with the empirical results. As our goal was to find the simplest model system that agreed with experiment, we pursued further simplification by including only Fe2+ and the ligand (without the carboxylic acid; see above). The total energy (Fe2+ complexation and strain) calculated with this simpler model system still favored the pyrazole motif. This became the model system for subsequent investigations and will be the basis of the results reported herein.

Figure shows the results of the calculations on both benzimidazole analogs (1 and 2). The pyrazole motif (i.e., 2) has stronger Fe2+ complexation (by 1.78 kcal/mol) and a lower strain penalty (by 1.41 kcal/mol), which results in it being favored overall by 3.19 kcal/mol. The higher strain for the amide motif of 1 can be understood in terms of the bound-state conformation needing to be nearly planar for the portions of the molecule (e.g., benzimidazole and amide) that complex Fe2+, which results in an unfavorable dipole–dipole interaction between the NH groups of benzimidazole and the amide bond. Indeed, a similar disfavored interaction occurs for the pyrazole motif. However, it is to a lesser extent since the dipole–dipole interaction is now between the NH of benzimidazole and the CH of pyrazole.

Figure 2

Calculation scheme for 1 and 2. ΔΔEcomplex is the difference in Fe2+ complexation energies, while ΔΔEstrain is the difference in strain energies between the two analogs. In both cases, the energy differences favor the pyrazole motif for a total of 3.19 kcal/mol. The key portion of the tail motif is highlighted in blue; the remaining acid portion of the tail motif is highlighted in red and (as noted in the text) was not explicitly included as part of the model system. The arrows indicate the strain between adjacent regions.

Having satisfactorily reconciled the nature of the potency difference between compounds 1 and 2, we then focused the structure–activity relationship (SAR) development around the more potent benzimidazole-2-pyrazole carboxylic acid series.

Fe2+ Complexation Energy Calculations

Three separate minimizations are performed to obtain ΔEcomplex: (i) minimization of the Fe2+–molecule complex, (ii) minimization of Fe2+, and (iii) minimization of the (fully unconstrained) molecule. To obtain the ΔEcomplex (of formation) for a given complex, the final energies of Fe2+ and the molecule are subtracted from the final energy of the Fe2+–molecule complex. To obtain the ΔΔEcomplex of two different complexes, the ΔEcomplex of the lower energy complex is subtracted from the ΔEcomplex of the higher energy complex. By this definition, the ΔΔEcomplex is a positive value that indicates how much lower in energy the lower energy complex is versus the higher energy complex. All minimizations were performed using Jaguar (under Maestro Version 10.7.015, release 2016-3) with density functional theory (B3LYP) and the LACV3P* basis set. Further, symmetry was turned off, the grid density was set to “fine”, and the accuracy level was set to “accurate”. We note that to find the lowest energy state for a given species, it may be necessary to slightly vary the input conformation/geometry and to perform multiple minimizations. (For example, even for the Fe2+–molecule complex minimizations, we would vary the position of the molecule relative to Fe2+ and perform multiple minimizations.) We then deemed the final energy as the lowest energy value that converged.

Strain Energy Calculations

The molecule conformation obtained from the minimization of the Fe2+–molecule complex is extracted. The dihedral angle of the molecule is constrained to its original value, and a subsequent minimization is performed on the molecule alone. A fully unconstrained minimization of the molecule is also performed. To obtain the ΔEstrain of a given molecule, the final energy of the unconstrained (lower energy) conformation is subtracted from the final energy of the constrained (higher energy) conformation. By this definition, the sign of ΔEstrain is positive. To obtain the ΔΔEstrain of two different molecules, the ΔEstrain of the lower strained molecule is subtracted from the ΔEstrain of the higher strained molecule. By this definition, the ΔΔEstrain is a positive value that indicates how much lower in strain the lower strained molecule is versus the higher strained molecule. All minimizations were performed using Jaguar (under Maestro Version 10.7.015, release 2016-3) with density functional theory (B3LYP) and the LACV3P* basis set. Symmetry was turned off, the grid density was set to “fine”, and the accuracy level was set to “accurate”. We note that to find the lowest energy state for a given molecule, it may be necessary to slightly vary the input conformation and to perform multiple minimizations. (Indeed, even for the dihedral-constrained molecule minimizations, we would slightly displace the dihedral from its starting value, upon which it obtained the desired constrained value upon completion of the minimization.) We then deemed the final energy as the lowest energy value that converged.

Substitution Effects on Fe2+ Complexation and the Preferred Tautomer

The effect of 5,6-disubstitution of the benzimidazole ring on the Fe2+ complexation energy for the benzimidazole-2-pyrazole carboxylic acid series was investigated. For all analogs of interest, the Fe2+ complexation for each of the two orientations was calculated.

Table shows the lower energy orientation (with respect to Fe2+) for several 5,6-disubstituted analogs along with the energy difference between the two orientations (the negative sign indicates that the orientation shown is preferred according to the model). An interesting trend emerges: the more electronegative substituent of the two faces Fe2+. Moreover, these calculations reveal which benzimidazole nitrogen interacts with Fe2+ and, in turn, which nitrogen bonds with hydrogen; as such, we have deemed the resulting energy difference the preferred tautomer energy (ΔEtautomer).

Table 1

Experimental pIC50 Values and Calculated Values for ΔEtautomer and pKa (DMSO) for Several 5,6-Disubstituted Benzimidazole-2-pyrazole Carboxylic Acid Analogs (2a–2f)a

compound	R1	R2	ΔEtautomer (kcal/mol)	pKa (calcd in DMSO)	pIC50 (exp)
2a	Cl	pyrrolidine	–2.74	12.9	6.65
2b	Cl	MeO	–2.35	12.6	6.50
2c	F	Cl	0.65	11.0	6.70
2d	CF3	Cl	–2.08	10.2	7.14
2e	CN	Cl	–1.24	9.7	7.40
2f	NO2	Cl	–1.74	9.1	7.50

The expected (from the ΔEtautomer) binding orientation with respect to Fe2+ is shown. Highlighted in blue is the pKa proton. The acid group is highlighted in red to remind the reader that this motif was removed to further simplify all calculations.

Compared to docking alone, this calculation provided important information that allowed the determination of binding poses with greater certainty. For example, in Figure , we see the crystal structure of 2c, which is in excellent agreement with the binding pose that was predicted (prior to obtaining the crystal structure) using the combination of docking and Fe2+ complexation energy calculations. In general, with respect to its interaction with Fe2+ and other key interactions (shown as dashed lines), this binding pose is representative of the benzimidazole-2-pyrazole carboxylic acid series. Further, as seen in Figure , these are the same interactions made by 1.

The Water Model

Although our crystal structures revealed a conserved cascade of water in the active site, we chose to focus on the single water molecule (of this cascade) closest to the protonated nitrogen of benzimidazole since (as noted earlier) it showed a bridging water interaction to Y303. The first step in constructing our water model was to determine which nitrogen of disubstituted benzimidazole was protonated; in other words, we first calculated the preferred tautomer as previously described.

Thereafter, Jaguar[10−12] was used to calculate the pKa of this key proton. As part of the Fe2+ complexation energy calculations, it was necessary to optimize the Fe2+–ligand complex. From this calculation, we took the ligand conformation alone and used it as the starting conformation in Jaguar. (We note that Jaguar does not currently allow one to fix the input conformation and the inclusion of Fe2+ as part of the pKa calculation is not recommended in general.[13])

Table shows the calculated pKa (in DMSO) versus the experimental pIC50 values. In general, as the benzimidazole proton becomes more acidic (as the pKa decreases), the pIC50 increases. We hypothesize that this proton is interacting with the oxygen in the water molecule to form a hydrogen bond. Thus, as the pKa decreases, it becomes a better donor, which results in potency improvement. In Figure , the correlation between the calculated pKa (DMSO) and pIC50 values is plotted. The good correlation (R2 = 0.870) provided a reliable way to easily obtain potency gains.

Figure 4

The correlation of the experimental pIC50 with the calculated pKa (DMSO) of the key proton of the benzimidazole series (R2 = 0.870).

The Cinnoline Series: A Prospective Model

As the project progressed, it became desirable to find other viable analogs and move beyond the benzimidazole-2-pyrazole carboxylic acid series. This also afforded the opportunity for the model to make prospective predictions.

From the literature,[14] 4-hydroxycinnoline-3-carboxamide was known to give good inhibition of prolyl hydroxylase activity, which led us to consider its viability for PHD2. Both the carboxamide and pyrazole carboxylic acid analogs (3 and 4, respectively) were interrogated in the modeling paradigm (docking, Fe2+ complexation energy, ΔEtautomer, and strain).

Figure displays the results of the Fe2+ complexation energy and strain calculations on both 4-hydroxycinnolines motifs. Unlike the benzimidazoles, for the 4-hydroxycinnolines, it is the amide motif that has the stronger Fe2+ complexation (by 2.13 kcal/mol) and lower strain penalty (by 0.629 kcal/mol), which results in it being favored overall by 2.76 kcal/mol. The physical origins of the strain difference are similar as noted for the benzimidazoles, in which the portions of the molecule (e.g., cinnoline and pyrazole) that complex Fe2+ will be nearly planar. Whereas for the 4-hydroxycinnoline-3-amide motif, this results in a (favorable) intramolecular hydrogen bond interaction between the oxygen atom of the cinnoline hydroxyl and the NH of the amide; the pyrazole motif is unable to make such an interaction.

Figure 5

Calculation scheme for 3 and 4. ΔΔEcomplex is the difference in Fe2+ complexation energies, while ΔΔEstrain is the difference in strain energies between the two analogs. Unlike for the benzimidazoles, for the cinnolines, the energy differences favor the amide motif for a total of 2.76 kcal/mol. The key portion of the tail motif is highlighted in blue; the remaining acid portion of the tail motif is highlighted in red and (as noted in the text) was not included as part of the model system. The arrows indicate the strain between regions, whereas the dashed lines indicate a key hydrogen bond.

Both compounds were subsequently made, and as predicted by the model, the amide analog (3) was more potent (pIC50 = 6.9 versus 5.0 for the pyrazole analog (4)). That the pyrazole motif did not provide the potency gains seen for the benzimidazoles is a striking reminder of a non-additive SAR. In other words, whether the pyrazole or amide motif is preferred depends on what it is “paired” with and whether this resulting molecule can achieve the near-planar binding conformation necessary to complex Fe2+ without incurring too much strain. The ability of the model to predict this became a powerful tool in our further design efforts.

In Figure , we show the crystal structure of 3 in the active site of PHD2. Once again, the modeling paradigm correctly predicted the binding pose. The key interactions are like those seen before for the benzimidazole analogs, where it makes a salt bridge with R383 via the acid moiety and a bidentate interaction with Fe2+ via the N1 nitrogen atom of the cinnoline and the oxygen atom of the amide. However, unlike the benzimidazole series, which make a bridging water interaction to Y303 via the key water molecule, the 4-hydroxyl of the cinnoline displaces that water molecule completely and directly forms a hydrogen bond with Y303.

Figure 6

Crystal structure of 3 in the active site of PHD2 (PDB ID: 6NMQ). As before with 2c, the modeling paradigm correctly predicted the binding pose prior to obtaining the crystal structure.

Conclusions

The modeling approaches most often used in structure-based drug design are usually crystallographic analyses and docking. Indeed, these approaches are paramount, but moving beyond these qualitative approaches means systematically interrogating particular ligand interactions at a higher level of theory. Moreover, ligand binding strain is often a hidden culprit in ligand potency and should also be accounted for explicitly in computational models.

In the case of PHD2, devising computational models for the Fe2+ complexation and key water interaction (along with the ligand binding strain) allowed us to account for potency differences seen within both the benzimidazole (retrospectively) and 4-hydroxyl cinnoline (prospectively) series. Whereas the pyrazole motif was favored over the amide in the benzimidazole series for its potency gains, the reverse scenario was true in the 4-hydroxycinnoline series. This was a stark reminder of the often non-additivity effect that can frustrate SAR interpretation between different chemical series. The model’s ability to differentiate this effect among diverse series provided substantial guidance in future design efforts.