Abstract

The effect of weak intermolecular interactions on the binding affinity between ligand-protein complexes plays an important role in stabilizing a ligand at the interface of a protein structure. In this review article, we will explore the different ways of taking into account these interactions, mainly intramolecular hydrogen bonds, in docking calculations. Their possible limitations and their suitable application domains are highlighted. Inspection of the outliers of this study probed very stimulating, as it provides opportunities and inspiration to medicinal chemists, being a reminder of the impact that minimal chemical modifications can have on biological activities.

1. Introduction

Non-covalent interactions, especially hydrogen bonds (HBs), are known to be responsible for the conformation and 3D structure of biomolecules. Existence of protein’s secondary structure ¹, DNA and RNA are due to this kind of force. Other non-covalent interactions as halogen bonds ^{2, 3, 4}, pnicogen bonds ^{5, 6}, and tetrel bonds ^{7, 8} can contribute as well to the stability of certain molecular conformations.

Hydrogen (H)-bonds are of immense importance in biological systems ^{9, 10, 11} and ubiquitous in nature, playing an important role in protein folding ¹², protein-ligand interactions ¹³, and catalysis ^{14, 15} Despite extensive investigations, there remain many challenges that prevent us from completely understanding how H-bonds modulate molecular function. For that reason, hydrogen bonding has been one of the most interesting areas of research, primarily due to its crucial role in governing the shapes, properties, and functions of biomolecules. It affects many chemical, physical, and biological systems and processes ^{16, 17, 18}. Thus, methods that allow the prediction of the hydrogen bond strength directly from the molecular structure are particularly helpful in understanding and anticipating the molecular disposition for anticipated interactions.

In biological systems, an H-bond competing process is always present with water. Because bulk water interferes with reversible biological processes and enthalpy-entropy compensation occurs during H-bond formation, the mechanisms and the extent to which H-bonds contribute to molecular function are not well understood. In particular, whether H-bonds regulate receptor-ligand binding remains a long-standing problem with poorly defined mechanisms ^{19, 20, 21, 22, 23}.

When studying non-covalent interactions in solution, desolvation becomes important as the interaction between substrate and solvent compete with the noncovalent interaction ²⁴. In order for an interaction to occur between two molecules in solution they must first be desolvated (Figure 1). If the solvent interacts strongly with the substrate the non-covalent interaction is destabilized. It is therefore important that π-π interactions are studied in a range of polar and non-polar solvents to fully understand them. However, due to the flat π-electron surfaces of aromatic molecules, solvophobic forces favor π-stacking geometries which allow for maximum overlap and these are not commonly observed for π-π interactions. As a consequence solvophobic forces do not determine the geometry of these interactions.

The importance of the hydrogen bond to drug design is well recognized. Hydrogen bonds are not only crucial in dictating the orientation of an inhibitor binding in the receptor but also contribute importantly to binding affinity. Hydrogen bond capacity is an essential factor in the strategy of bioisosterism for drug design and optimization. When a bioisostere is used to replace an existing moiety of a compound, the replacement has to match the hydrogen-bonding characteristics of the parent and would preferably further improve upon compound properties including binding potency. There is strong evidence indicating that the strength of different hydrogen bond donors and/or acceptors varies significantly. The hydrogen-bonding constants of commonly encountered hydrogen-bonding groups have been measured to vary over more than 3 orders of magnitudes ²⁵. Furthermore, in the drug optimization process, it is often observed that electron-withdrawing or donating substituents have opposite effects on the activity of the compounds. Some of the underlying causes of such structure-activity relationships (SAR) could be traced to a modulation of the hydrogen-bonding interaction of the inhibitor with the receptor ²⁶.

H-bonds are generally considered to be facilitators of protein-ligand binding ^{13, 27}. However, introducing H-bond donors or acceptors to establish stronger protein-ligand interactions often results in the absence of net gain in binding affinity ^{19, 28}. Rather than targeting protein-ligand interactions per se, H-bonds are also reported to promote ligand binding affinity by displacing protein-bound water molecules into the bulk solvent ^{29, 30, 31, 32}. Contrasting H-bonding mechanisms are also evident in enzyme catalysis where the effects of H-bonds on the free energy barrier reduction of enzymatic reaction are identical to their role in protein-ligand binding. Whether electrostatic (H-bond) interactions represent the major origin of enzymatic catalytic power is still under debate ^{33, 34, 35}. We recently reported that accurate quantification of the free energy contribution of H-bonds to both enzymatic reactions and the corresponding reference reactions in aqueous solution is vital for exploring the origin of enzymatic catalytic power ³⁶. A deeper understanding of the effects of water H-bond interactions on biological processes is therefore needed to advance our appreciation of how such systems are regulated and to facilitate lead compound design without compromising ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties or escalating drug development costs ³⁷.

Clearly, the variations in hydrogen bond strength could be utilized in drug design. A fundamental understanding of hydrogen bond interactions will greatly facilitate such efforts. In a comprehensive monograph on the subject ¹⁶, Jeffrey and Saenger summarized the known characteristics of hydrogen bonds, some of which are particularly instructive to a theoretical investigation: (1) hydrogen bonds are not properties of atom pairs but are dependent on the pair of atom groups that forms the extensive donor and acceptor subunits. ^{16, 25}; (2) the major component of hydrogen-bonding interaction is electrostatic ^{16, 38, 39}; and (3) hydrogen bonds are soft interactions, and hydrogen bond lengths and angles fluctuate according to local environments ^{16, 40, 41}.

As H-bonds play an important role for the ligand-protein binding process, being responsible for direction and recognition of substrates and modifying the affinity to their binding partners ⁴², it is important to take into account the fact that the development of new drugs is a very time-consuming process requiring a huge financial investment ⁴³ and it is aiming for a high affinity and selectivity of the ligand’s binding to its target ⁴⁴. To achieve high drug binding affinities either binding enthalpy and binding entropy must contribute in a favorable way to binding since from a thermodynamic point of view the binding constant K_a is defined as

(1)

where T is the absolute temperature, R is the gas constant, and ΔG is the change in Gibbs free energy that in its turn is given by the equation

(2)

Here ΔH and ΔS are the changes in binding enthalpy and binding entropy, respectively. Related to formula 2 it is clear that strong binding affinities can be reached either by a more negative ΔH, a more positive ΔS, or a combination of both ⁴⁵. At physiological conditions a change of 1.4 kcal/mol in Gibbs free energy ΔG leads to a 10-fold higher or lower equilibrium constant Ka as a result of the exponential relationship between these two thermodynamical properties ⁴⁶ (see eq. 1). The binding enthalpy itself depends on the interaction forces (van der Waals, hydrogen bonds, etc.) between the ligand and the target protein, while the binding entropy is made up of two main contributions: changes in solvation and conformational entropy. This is evidence of the fact that drug affinity is related not only to structural behavior (enthalpy) but also the dynamics (entropy) of the interacting species ⁴⁷. The simultaneous optimization of the enthalpy and entropy term is though a difficult goal to achieve, because it implies the overcoming of the so-called ‘enthalpy/entropy compensation’ phenomenon. This effect consists in an entropy loss contemporary to an enthalpy gain that can be understood when a ligand binds to its protein by establishing favorable interactions like hydrogen bonds but losing conformational freedom. Further the optimization of the binding enthalpy depending on various forces is very difficult, while the same optimization process for the binding entropy is easier to obtain because of its primary dependence on the hydrophobic effect. One strategy in drug design is therefore the generation of hydrophobic and conformationally constrained ligands ⁴⁸.

Their binding affinity is entropically dominated, but binding enthalpy often shows an unfavorable contribution. Those ligands being constrained in their conformation cannot easily respond to binding site geometry changes, and so they are highly susceptible to drug resistance mutations or genetic polymorphism naturally occurring ⁴⁸. Also such compounds with entropically dominated binding show a significant improvement in binding affinity when an unfavorable binding enthalpy is eliminated. The importance of optimizing binding enthalpy during drug design is consequently clear. In this context hydrogen bonds play a key role in the gain of enthalpy. They are also crucial to improve selectivity as they are determined by strict geometric and distance constraints ⁴⁹. An aim in drug design is therefore to modify the hydrogen bonds strength between a ligand and a protein in order to achieve higher binding affinities but without any negative influence on other interactions contributing positively to binding established ⁵⁰. One way to accomplish this intention is to alter the chemical properties in the neighborhood of an existing H-bond by adding functional groups or to modify existing groups ^{51, 52}.

At the same time, the study of intramolecular interactions is very important in the design of pharmaceutical drugs, particularly in the context of conformationally flexible molecules. Conformation-controlling intramolecular interactions in drug molecules have a direct influence on the binding modes of the drugs with the respective targets ^{53, 54, 55}. In particular, intramolecular peri-interactions have been widely studied in the literature in naphthalene and other related systems ^{56, 57, 58, 59, 60}. For that reason, very recently Caron et al. stated the necessity to implement intramolecular hydrogen-bonding considerations in drug discovery ⁶¹.

2. Nature of Hydrogen Bonds

The hydrogen bonds have some interesting features like: (i) low enthalpy formation, about 20% of the chemical bond enthalpy; (ii) total unspecificity; (iii) the association, or dissociation, of molecules through HB are fast enough to permit to check their formation and to correct misformations. The large number of HB, that are generally present between water and biological materials, result in the high specificity of the HB networks which include intra and intermolecular bonds.

Conformational preferences can cause non-contiguous atoms within an isolated molecule to become similarly close neighbors. These spatial arrangements may be driven by favorable electrostatic interactions or by the special case where three of such atoms form a so-called “hydrogen bond” (H-bond). Although the situation becomes more complicated when the molecular structure is considered within a solution environment, these same two factors remain important to also drive additional intermolecular interactions now possible between solute molecules themselves and with the solvent molecules as partners. Focusing herein on hydrogen bonding, it can be noted that, despite a decades-long endeavor to define the H-bond, this key arrangement still cannot be considered to be resolved with full consensus. The 2011 IUPAC recommendations provide a definition ⁶² that can be used as the basis for critical evaluation.

The recommendations state: “The hydrogen bond is an attractive interaction between a hydrogen atom from a molecule or a molecular fragment X–H in which X is more electronegative than H, and an atom or a group of atoms in the same or a different molecule, in which there is evidence of bond formation. A typical hydrogen bond may be depicted as X–H…Y–Z, where the three dots denote the bond. X–H represents the hydrogen bond donor. The acceptor may be an atom or an anion Y, or a fragment or a molecule Y–Z, where Y is bonded to Z. In some cases, X and Y are the same. In more specific cases, X and Y are the same and X–H and Y–H distances are the same as well leading to symmetric hydrogen bonds. In any event, the acceptor is an electron rich region such as, but not limited to, a lone pair of Y or π-bonded pair of Y–Z. The evidence for hydrogen bond formation may be experimental or theoretical, or ideally, a combination of both. Some criteria useful as evidence and some typical characteristics for hydrogen bonding, are not necessarily exclusive. The greater the number of criteria satisfied, the more reliable is the characterization as a hydrogen bond”.

The first criterion for a hydrogen bond claims: “The forces involved in the formation of a hydrogen bond include those of an electrostatic origin, those arising from charge transfer between the donor and acceptor leading to partial covalent bond formation between H and Y, and those originating from dispersion.” It reveals from the specification in footnote F2 that “Attractive interactions arise from electrostatic forces between permanent multipoles, inductive forces between permanent and induced multipoles, and London dispersion forces. If an interaction is primarily due to dispersion forces, then it would not be characterized as a hydrogen bond.”

In a special case, Y symbolizes an aromatic ring with its electron cloud favorably interacting with a positively polarized H atom. This sort of hydrogen bond is called an H…π interaction. For X–H…Y with X = C or with X, Y = S or P, as well as for the H…π interaction, the role of the dispersion forces increases in comparison to the cases where the H-bond formation is principally related to electrostatic and charge-transfer effects.

The third criterion (E3) on the list of the IUPAC recommendations says: “The X–H…Y angle is usually linear (180°) and the closer the angle is to 180°, the stronger is the hydrogen bond and the shorter is the H…Y distance.” Two important footnotes were added to this criterion. “The X–H…Y hydrogen bond angle tends toward 180° and should preferably be above 110° (F4).” “Historically, the X to Y distance was found to be less than the sum of the van der Waals radii of X and Y, and this shortening of the distance was taken as an infallible indicator of hydrogen bonding. However, this empirical observation is true only for strong hydrogen bonds. This criterion is not recommended. In most cases, the distance between H and Y are found to be less than the sum of their van der Waals radii. It should be noted that the experimental distances are vibrational averages and would differ from such distances calculated from potential energy minimization. (F5)”.

Thus, as revealed by the quoted text, no H…Y or X…Y distance has been strictly defined for the distance of a H-bond, nor has a strict lower limit for the X–H…Y angle has. On the other hand, the X…Y distances for the different intramolecular H-bonds could represent borderline cases with values equal or slightly larger than the sums of the van der Waals radii. Likewise, in cases when a H-bond can form a five-member ring arrangement (Figure 1), the X–H…Y bond angles could be close to or even less than 110.

In aqueous solutions, the O (solute)…O (water) and N (solute)…O (water) radial distribution functions show their first minima at up to 350 pm ^{63, 64}. This value has been accepted as the boundary of the first hydration shell around the polar sites of solutes. This, however, does not mean that intermolecular H-bonds would be expected with X (solute)…O (water) separation up to 350 pm in solution. Analyses of Monte Carlo results (see below) always point out that the number of the solvent molecules engaged in H-bond(s) to the solute is smaller than the total number of the solvent molecules in the first hydration shell(s) around the polar site(s).

Thousands of different proteins in the cells, precisely codified by the genes, realize specific functions. The individuality of the proteins is directly related with their three-dimensional structure that provides the ideal conditions for realizing correctly their functions ⁶⁵. However, in spite of the importance of the proteins in the cellular life, the role of the supporting medium, in particular of the solvent molecules, i.e. the water molecules, is equally important because they constitute the most part of the cell, about 70%, and because the contribution of the water molecules is essential in the peptide activities ⁶⁶. The water molecules ^{67, 68}, beyond the small mass, present high multipolar moments that contribute to the formation of hydrogen bonds, HB. Consequently, the water molecules perform an essential structural role in the organization and activity of the biological medium. Many of the protein activities depend on the protein stability, on associations with other proteins or ligands while the catalytic activity depend on the structure, on thermodynamic and dynamic properties, properties that are deeply influenced by the solvent. All the most important biological molecules like peptides, saccharides, nucleic acids have the common feature that they contain hydrogen-bonding functional groups ⁶⁵.

The hydrogen bond (HB) is an attractive interaction between a hydrogen bond donor (HBD), XH, and a hydrogen bond acceptor (HBA), Y, in the same (intramolecular HB) or in a different (intermolecular HB) molecule. It is usually represented as XH···Y, where X is an electronegative atom such as O, N or F, and Y can be any electron-rich atom ⁶⁹. Although Y atoms bearing lone pairs of electrons provide the most favorable site for HB formation, hydrogen bonding through π-electron systems can occur. The HBD and HBA groups that are commonly seen in biomolecules ¹⁶ are shown in Table 1.

Strictly speaking, the HB strength must be measured by the binding energy. Most HBs have a binding energy (1 cal = 4.184 J) in the range 2–20 kcal mol^{−1 70}. For example, the HB in the water dimer has a binding energy of 5.4 ± 0.7 kcal mol^{−1 71}. It is not easy to measure binding energies accurately and very few results are available for hydrogen bonded complexes.

Most chemists and biochemists characterize, instead, the HB strength from the changes in thermodynamic state functions pertaining to the formation of the hydrogen bonded complex XH· · ·Y from species XH and Y in solution ⁷². The Gibbs energy change ΔG of Reaction 3 is related to the formation or equilibrium constant K, the formation enthalpy ΔH and the formation entropy ΔS of the complex by Equation 4.

(3)

(4)

Many physical properties familiar to the medicinal chemist, such as the hydrophobic, electronic and steric substituent constants of the Hansch equation ⁷³, as well as the biological activities themselves, are related to the Gibbs energy. It seems logical, therefore, to set up scales of HB strengths using Gibbs energies rather than enthalpies.

Following the pioneering work of Taft and co-workers ⁷⁴, we have constructed a quantitative scale of HBA strength on these lines. We defined a quantity pK_HB as log K for the 1:1 complexation of a reference HBD, 4-fluorophenol, with HBAs, in carbon tetrachloride at 298 K (Equations 5–8). These pK_HB values, or the related ΔG_HB values (expressed on the molar concentration scale), then represent the relative HBA strength of molecules Y towards the reference HBD 4-fluorophenol.

(5)

(6)

(7)

(8)

It is important to note that pK_HB is defined as the antilog of the dissociation equilibrium constant so that it is in accordance with the definition of the pK_a scale. The most basic and the strongest hydrogen bond acceptor compounds have respectively the largest pK_a or pK_HB values.

A hydrogen bond between XH and Y consists in the sharing of a proton between Y and the anion X⁻ derived from XH, and therefore constitutes a step before proton transfer ⁷⁷:

(9)

Therefore, HBAs are bases, HBDs are acids and the pK_HB scale is a basicity scale. We shall show later that the proton sharing basicity (pK_HB) differs significantly from the proton transfer basicity (pK_a). Thus, the bases dimethylacetamide, dimethylsulfoxide and hexamethylphosphoramide are much weaker proton acceptors but significantly stronger hydrogen bond acceptors than pyridine or trimethylamine. This is particularly important because both proton transfer basicity and hydrogen bond basicity play important roles in biological systems, and it is often necessary to distinguish between the phenomenology attributable to both types of interactions.

The strength of HBs depends on the HBA and HBD strengths. Since HBDs of biological interest are mostly confined to OH and NH groups and since HBAs vary much more in character, we have first covered the scale of HBA strength. For HBD strength scale, we presently refer the reader to the log K^A_H ⁷⁸ or log K_α scales ²⁵. The validity range of the pK_HB scale, constructed towards a phenol in a medium of low reaction field, extends over many OH and NH donors and to higher reaction fields than the definition one. Thus, the pK_HB scale gives the medicinal chemist an extensive database for QSAR studies and, via structure–basicity relationships ^{79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101}, structural tools for optimizing the HBA strength of drugs. The pK_HB scale can also support spectroscopic studies of drug interactions and computer modelization of drugs through its correlations, generally family-dependent, with a number of spectroscopic, geometrical, electrostatic and thermodynamic parameters of the hydrogen bond, namely NMR ⁷⁴, UV and IR shifts ^{79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101} upon hydrogen bonding, hydrogen bond length ¹⁰², electrostatic potential ¹⁰² near the HBA lone pair, and enthalpy of HB formation ^{72, 102}.

In intramolecular hydrogen bonds, the geometry for both the H…Y distance and the X–H…Y angle is primarily determined by the covalent structure of the molecule. While three-atom hydrogen-bonded rings are extremely rare, the four-atom substructures (e.g., carboxylic group, amides) deserve special consideration. In most cases, a H-bond can be expected if the system can form a five to seven-member ring, including arrangements utilizing a polar H. Prototypes are indicated in Figure 2 and typical representatives of five- and six-member rings are shown Figure 3 and Figure 4 ¹⁰³. Seven-member rings can be formed for γ-substituted carboxylic acids, 1,4-disubstituted butanes with OH and/or NH₂ substituents. Larger rings are probably not stable.

When a molecule dissolves, a close molecular environment is encountered that is in contrast to the most frequently applied ideal-gas model, where no potential energy interaction is considered even through the collisions of the molecules. Although the solute-solvent interactions are substantial, the effect of a non-polar or only slightly polar solvent (CCl₄, CHCl₃) on the molecular geometry is generally small ^{104, 105}. The geometric effect could be, however, large when a solute with an intramolecular H-bond in the gas phase dissolves in a protic solvent such as water or methanol, which has both proton donor and acceptor sites. In this case, the X–H…Y intramolecular H-bond may collapse while solute-solvent H-bonds are formed using the free XH and Y sites.

The weakest point of the continuum dielectric solvent model is that the above solute-solvent H-bond(s) are only implicitly mimicked by polarization of the solvent and concomitant appearance of surface charges on the inner surface of the cavity: Negative surface charges opposite to a polar hydrogen and positive ones in the lone-pair regions of the solute’s oxygens and nitrogens. Although this response is qualitatively correct, the calculated solute-solvent stabilization energy is underestimated ^{106, 107}.

Thus, for proper calculation of the free energy changes when a polar solute with or without internal H-bond(s) dissolves in a protic solvent, explicit consideration of the solute-solvent intermolecular H-bonds becomes necessary.

This requirement can be largely satisfied by adopting the supermolecule + continuum approach, where the solute is surrounded by a number of explicit solvent molecules. The solute and the explicit solvent molecules mimic the H-bonds in the first solvation shell within the cavity carved in the continuum solvent. The critical question then becomes, how many explicit solvent molecules are to be considered.

For constructing the starting geometry of a supermolecule, knowledge of microsolvated solute structures is very helpful. In these systems, the central, polar molecule with or without an intramolecular H-bond is solvated by a few solvent molecules. Locations of the solvent molecules (water, methanol) indicate the most preferable solvation sites of the solute with a hydrogen donor/acceptor solvent.

Useful information can be obtained from experimental gas-phase hydration/solvation studies augmented with theoretical calculations ^{108, 109, 110, 111, 112, 113, 114, 115, 116} or specific theoretical calculations for hydrated amino acid side chains, nucleotid base and sugar models ^{117, 118, 119, 120, 121}.

There are many definitions of the nature, diversity, and range of the strength of hydrogen bonds ^{62, 122, 123, 124, 125} resulting from the covalent three-center four-electron bond and the electrostatic attraction and repulsion interactions between electron-rich groups. Intramolecular and intermolecular hydrogen bonds can also function as cooperative HB, three-center donor (bifurcated) HB and three-center acceptor (anticooperative) hydrogen bonds ^{17, 122, 126, 127, 128}. A special class of intramolecular hydrogen bonds predicted to be very strong because the neutral donor and acceptor atoms are linked by a system of π-conjugated double bonds was defined by Gilli ^{129, 130} as a resonance-assisted hydrogen bond (RAHB). The investigations by many authors ^{125, 131, 132} have established that RAHB systems by their aromaticity numerical descriptors (aromaticity indices ¹²⁵) fulfill (at least partially) the conditions of the aromatic π-electron delocalization; therefore, they can be treated as quasi-aromatic.

One of the first observations of hydrogen bonding noted that the O···O distance is, with certain exceptions, a reliable indicator of intramolecular HB strength ^{133, 134, 135, 136}. The synthetic parameter Q introduced by Gilli and co-workers ^{129, 130, 137} for the description of the RAHB intramolecular H bonds in the enols of β-diketone (π-delocalization of the O=C−C=C−OH enolone fragment) can be calculated with considerable accuracy by considering the C−O, C−C, C=C, and C=O crystallographic distances in the HB rings ^{138, 139}. However, correlations between the Q parameter and the H bond energy ¹⁴⁰ are not always present; the same is true for the various modified Lippincott−Schroeder models ^{126, 133}.

Hydrogen bonding may be considered a special class of Lewis acid−base interactions ^{141, 142}. Many examples define the role of the hydrogen bond based on the acidity−basicity concept in the properties of interest for drug design. Laurence ¹⁴¹ described in detail the hydrogen bond basicities of functional groups relevant to medicinal chemistry. Recently, the empirical appreciation of the bond strength according to the pK_a(H₂O) equalization principle was used for the prediction of H bond strengths from the acid−base parameters of the interacting partners ^{142, 143}.

3. Types of Hydrogen Bonds

Looking first at hydrogen bonding as a whole, there are three fundamentally different types of hydrogen bonds that each employ the same weakly acidic X–H proton donor groups (X=N, O..) but that differ in the nature of the proton acceptor. In the classical hydrogen bond (X–H···l.p type; l.p = lone pair), the proton acceptor is a lone pair of a weakly basic, electron rich element, typically N, O or halide ion. This type has been recognized since the 1930s and is of most importance in biochemistry ⁶⁵. Very recently, p-electrons, such as those of arene rings or CC multiple bonds, have been shown ^{144, 145} to be able to act as weak proton acceptors in hydrogen bonding. This X–H···π type is weaker than the X–H···l.p type, probably because π bonding electrons are in general much less basic than lone pairs.

In view of the decrease in hydrogen bond strength on going from X–H···l.p to X–H···π types, one would have expected to find that if any X–H···σ type of hydrogen bond existed it would be even weaker than the X–H···π type and so neither be readily detectable nor have significant effects on physical properties, σ-bonding electrons being even less basic than π-bonding ones. In fact, it is now clear ^{146, 147, 148, 149} that such hydrogen bonds can in certain cases be much stronger than the X–H···π type and comparable in strength with the classical X–H···l.p type. This situation can be found in metal complexes, so they are not included in the scope of this review.

Nishio et al. first introduced the idea of aromatic rings being involved in a weak form of hydrogen bonding ^{150, 151, 152, 153, 154, 155}. They highlighted how C-H/π hydrogen bonding can account for the folding tendency in a wide range of small organic molecules and how it plays an important role in the conformational behavior of organic molecules. The enthalpy for a single unit of a C-H/π interaction is small, approximately 0.5 to 2.5 kJ mol-1, when an alkyl or aromatic C-H is involved ¹²². However the total enthalpy can become sizeable when several C-H groups can simultaneously participate in interactions with π groups. As multiple C-H/π interactions formed between side chains can be significant, they may be considered as one of the driving forces to constrain a peptide conformation and consequently direct specific conformation in many proteins ¹⁵⁶. Using the crystal structure database Nishio et al. investigated the C-H/π interactions in peptides ^{157, 158}. Evidence to indicate that 42% of the structures studied exhibited such aromatic interactions was found. As a consequence the importance of this interaction for the folded conformations of peptides was realized.

In 1993, Sakai reported a computational study of the binding energies and structures of benzene-methane complexes and suggested that the dispersion force was the most important for C-H/π hydrogen bonding (Figure 5) ¹⁵⁹. For the most stable methane-benzene complex 11, the orientation of the methane is above the plane with one hydrogen directed towards the center of the aromatic ring with an angle approaching 180°, thus although unimportant when considering the binding energy, electrostatics determine the geometry of C-H/π complexes.

However, in contrast to conventional hydrogen bonds, the directionality of C-H/π bonds is weak. The binding energy of the most stable complex (11) was calculated to be -0.57 kcal per mol. More recent theoretical studies agree with these findings ^{160, 161, 162, 163, 164}. However, this only holds for a typical sp³ C-H/π hydrogen bond. For sp hybridized CH groups the contribution from electrostatic energy becomes significant, and as a consequence such interactions are far more similar in nature to conventional hydrogen bonds than sp³ C-H hydrogen bonds ¹⁶⁵. In essence, the proportion of electrostatic forces is dependent on the hybridization of the carbon atom in the C-H bond. Computational studies have also been used to determine the interaction energies and geometries of the benzene-acetylene and benzene-ethene complexes. For the acetylene-benzene complex the most stable geometry 20 is a T-shaped stacking between the acetylene and benzene, with a distance of 3.5 Å between the carbon atom of the acetylene and the center of the ring. In contrast, the ethene-benzene complex has two stable geometries 21 and 22, and the more stable is complex 22.

In order to probe the hydrogen bond character of C-H/π interactions, Nishio et al. studied the substituent effect on a range of aromatic molecules capable of forming intramolecular C-H/π interactions, by NOE enhancement ¹⁶². The effect of substituents is a useful probe for the hydrogen bond character of C-H/π interactions, as in the formate ester 23 (Figure 7) since, if it behaves like a conventional hydrogen bond, an electron donating substituent on the hydrogen bond acceptor should increase the interaction. From peptide studies it was found that C-H/π interactions are favored in 5- and 6-membered rings, and this prompted a further study of a series of organic molecules capable of forming 5-, 6- and 7-membered ring C-H/π interactions (Figure 8) ^{166, 157, 168}. Electron donating substituents on the aromatic ring should raise the energy of the highest occupied π-orbital, and electron withdrawing substituents on the carbon atom of the C-H donor should lower the C-H antibonding orbital and hence the energy gap of the interacting orbitals required for a C-H/π hydrogen bond, thus making the interaction favorable. It was proved that this combination of substituents leads to the largest ratio of the folded conformers.

In this form of non-covalent interaction, the π cloud of the aromatic system acts as an acceptor and X-H, in which X is either an oxygen, nitrogen or sulphur atom, acts as a donor ¹⁶⁹. Like the C-H/π interaction, this interaction can be considered as a “non-conventional” form of hydrogen bonding. Along with conventional hydrogen bonding, this non-conventional form has also contributed to stabilizing effects in structural biology ^{170, 171}. Early evidence for X-H/π hydrogen bonding comes from infrared spectroscopy studies conducted in the 1950s and 1960s by David et al. ¹⁷², Josien et al. ¹⁷³, Oki et al. ¹⁷⁴ and Yoshida et al. ¹⁷⁵. More recently, extensive searches of the Cambridge Crystallographic Structural Database have confirmed the presence of these interactions in a wide range of organic molecules ^{176, 177, 178}.

It was reported that the attraction between benzene and water or ammonia is weaker than for a conventional hydrogen bond. The interaction energy of the benzene-ammonia complex was found to be approximately 40% of the energy of a conventional hydrogen bond. As expected for hydrogen bonds, the order of magnitude of the E_tot values is in accordance with the electronegativity of the proton-donating atom, so the strongest interaction occurs in the benzene-water complex. The shortest distance between the center of the aromatic ring and the atom of the other molecule was also observed for the benzene-water complex. The results also show that in the case of the ammonia complex, the preferred geometry is with the ammonia above the center of the ring and the N-H bond in the perpendicular plane. This is indicated by the significantly larger interaction energy for complex 36, as compared to 39 and 40. For complex 36 a substantial attraction was found at intermolecular distances over 4.0 Å and the energy of the complex was still orientation dependent at these separations. This prompted them to suggest that long range and not short range interactions are the main source of attraction. These results show that the dominating forces are dispersion and electrostatics and not short range forces such as charge transfer.

4. Influence of Hydrogen Bonds in Binding

Hydrogen bonding is well recognized as one of the major noncovalent forces which play a prominent role in supramolecular and template chemistry and is a crucial issue in the study of biologically important molecules ^{16, 17, 18}. In aqueous solutions, the hydration of pendent H-bonding donor or acceptor groups retards the permeation of a drug molecule through biological membranes ^{179, 180, 182}. The predictable formation of certain intramolecular motifs has a significant influence on the ability of the molecule to engage in intermolecular hydrogen bonding ¹⁸³ and affects molecular properties including biological and pharmacological activities. Structural behavior of biochemical systems and the nature of hydrogen-bonded liquids can be better understood when the concept of single hydrogen bonds is adopted and their detection and quantification is possible. The energy of intermolecular interactions is simply measured as a difference between the energies of a H-bonded structure and its components. This approach cannot be used to study intramolecular interactions, and despite the fact that the hydrogen bonds are qualitatively well understood, ^{18, 122, 124, 125, 183} it is generally admitted that a simple method of calculation and quantitative data is needed ¹⁸⁴. The energy of such optimized structures stabilized by the intramolecular hydrogen bond was used to estimate its strength by different comparative methods: cis−trans analysis ^{132, 185, 186}, isodesmic reactions ¹⁸⁷ and conformational analysis with gradual or total rotation of the H-bond donor or acceptor ^{138, 188, 189}. The first of these methods is simplistic, because the energetic stabilization of the H-bonded conformer includes several contributions, such as the balance between attractive and repulsive terms, steric constraints of bulky groups, conjugation and other interactions, and this situation changes radically after turning out of the hydroxyl group. In the isodesmic method the subtle differences in the electron density of O−H···O=C-bonded molecules may not be noticeable. Consequently, the calculation of the energy of intramolecular hydrogen bonds EHB in malonal and its derivatives structures show a significant disparity of results for different methods from 2.77 to 43.63 kcal/mol ¹⁹⁰ and from 4.9 to 18.0 ¹⁹¹. Other procedures of the EHB calculations more recently described ¹⁹² may be applied only for some groups of compounds.

Therefore, different authors ^{135, 140, 186} calculating the energy using these methods called them as binding energy, energy increase upon H-bond removal or the energy difference between the closed and opened structures of the molecule instead of the strength of the H-bond.

Because H-bonding is one of the most important factors that determine binding of small molecules to target proteins. However, a paradox often arises. When compounds with nanomolar potency are found, they are usually very water insoluble. Such compounds have been termed high-affinity traps by Stella and Borchardt ¹⁸². Rules were developed by Lipinski ¹⁹³, the Rules of 5, to address the high-affinity trap paradox. These rules were gleaned from a retrospective analysis of marketed drugs, the assumption being that such drugs as an aggregate have the physicochemical properties that should serve as a baseline for molecules that are identified by HTSs. Two out of the four rules deal with H-bonding atoms: N and O. These rules place limits on these atoms and are general enough to address both permeability and the H-bonding cohesive component of high-affinity traps that produces strong lattice energy crystals. The rules themselves have changed very little over the past 15 years. This is a testament to their validity but also to the difficulty for humans to refine and adapt them to more specific situations. More specific rules regarding H-bonding have arisen out of the emerging area of crystal engineering.

Carbohydrates are biomolecules widely studied due to numerous applications, which range from pharmacological action, for example, the anticoagulant activity of heparin ¹⁹⁴, to technological interest in developing cellulose nanofibers ^{195, 196, 197} that can exhibit mechanical properties comparable to carbon nanotubes ¹⁹⁸. Such properties are related to the molecular architecture, which is governed by intra- and intermolecular interactions, such as hydrogen bond. In addition to carbohydrates, many other compounds and biological molecules exhibit hydrogen bonding as generally one of the main stabilizing interactions, generating great interest of the scientific community for new studies with emphasis on this type of interaction, for example, in supramolecular chemistry.

Druglike organic molecules almost invariably contain a number of functional groups capable of forming hydrogen bonds, rendering them soluble and giving them the ability to form specific interactions with their biomolecular targets. When a donor and an acceptor are in proximity on the same molecule, an equilibrium may exist between closed conformations in which an intramolecular hydrogen bond is formed, creating a temporary ring system, and open conformations in which the polar groups are exposed to solvent (Figure 10).

These sets of conformations are not only structurally distinct. It is intuitively clear that the closed forms, hiding polarity from the environment, should be more lipophilic and might display a higher membrane permeability, whereas the open forms should be more water-soluble. In drug discovery, it is therefore important to recognize the potential for intramolecular hydrogen bond formation and to be aware of its consequences.

Hydrogen bonding patterns are critical to the packing of proteins and thus to the formation of the cavities and crevices due to incomplete packing that are the binding sites for drugs. We now know quite a bit about the hydrogen bonding characteristics of drugs. The three dimensional characteristics of hydrogen bonding are well understood from analyses of hydrogen bonds in small molecule complexes in the Cambridge Structural Database and from hydrogen bonds in protein−ligand complexes in the Protein Data Bank. From an analysis of ligands containing multiple hydrogen bond donor groups it is apparent that within a single ligand it is very difficult to accommodate more than about two or three hydrogen bond donor groups with the precise geometry needed for the maximum enthalpic energy benefit arising from the optimal three dimensional positioning of a hydrogen bond ¹⁹⁹. This phenomenon is often discussed in the context of entropy–enthalpy compensation ²⁰⁰. As a result, as the number of potential hydrogen bond donors in a ligand increases it becomes increasingly likely that they will not contribute in a positive sense to ligand binding and likely will detract from ligand binding. This observation is very consistent with the common medicinal chemistry observation that it is difficult to improve potency by addition of hydrogen bond donor groups. The directionality of hydrogen bonding places restrictions on the number of hydrogen bonds in a ligand quite apart from the effect of hydrogen bond donor and acceptor groups on membrane permeability ²⁰¹.

H-bond interactions between a ligand and its target protein are known to play a critical role in determining the overall affinity between the two. Not surprisingly, such interactions also have been considered to significantly influence recognition of compounds as substrates by P-gp ^{202, 201, 202, 203, 204}. In many cases, especially when comparing different chemical series to each other, a count of the number of HBD and HBA groups or calculated TPSA may serve as simple tools to differentiate scaffolds for relative likelihood of P-gp recognition. However, a more rigorous approach also can be desirable for the description of molecules as has been advocated by Abraham et al. ²⁰⁵. Not all donors and acceptors are ‘equal’ in terms of the strength of the H-bond that they form with the corresponding partner and hence parameters that can differentiate or rank order relative H-bond strengths are valuable. Among some of the experimental sources of H-bond scales, measurement of the H-bond equilibrium constants of a diverse set of proton donors and acceptors by Morris et al. is considered to provide an important reference set ²⁵.

Although this type of hydrogen bonds are out of this review scope, it is worth to point out that intermolecular interactions play a crucial role in a variety of life and biomolecular processes, including drug bindings, enzymatic chemical reactions, immune responses, genetic transcription and translation, cell signaling transduction and protein folding ²⁰⁶. It is well established that the strengths of protein-ligand binding and protein solvent interaction are strongly dependent upon the detailed position and orientation of the interacting partners. Although modern X-Ray crystallization techniques readily provide us with three dimensional coordinates of biomolecules based on which various types of non-bonded interactions can be identified, simply summarize number of non-bonded interactions are inadequate for a mechanistic understanding of the binding event. Thus, it is necessary to perform high level quantum mechanical calculations to quantify the strength of intermolecular interactions and to provide mechanistic insights into binding, having into account the possibilities of some ligand groups to form both inter and intramolecular hydrogen bonds, depending on different specific circumstances.

Intramolecular hydrogen bonds (IMHBs) play an essential role in biochemistry and chemistry. They affect the electronic distribution, molecular geometry, shape, and conformation of systems as diverse as proteins, nucleic acids, catalysts, and materials. As such, IMHBs greatly impact molecular properties, function, and interactions ^{207, 208, 209, 210, 211}. For these reasons, the formation or disruption of IMHBs is used by medicinal chemists to modulate biological and chemical properties of interest. As an example, the bioactive conformation of a given ligand could be stabilized by IMHBs. This could reduce the conformational and translational entropy upon binding and result in stronger association ^{212, 213, 214, 215}. Furthermore, the accessibility of polar atoms in a molecule could be decreased if IMHBs are established. This may influence desolvation equilibria and facilitate the passage of molecules through low dielectric environments ^{216, 217, 218}. While geometric preferences for the formation of IMHBs have been described based on experimental evidence ⁶⁷, the presence of stabilizing IMHBs are normally debated based on thermodynamics grounds, especially in high dielectric environments. Furthermore, the effects originating from the formation or disruption of IMHBs have been difficult to predict from structure, increasing the uncertainty in robustly utilizing IMHBs as a design concept. The present study attempts to evaluate the impact of IMHBs on biological activity as a commonly used in vitro optimization variable.

While in general the bioactivity difference does not systematically and significantly deviate from ±0.5 log unit, it is interesting to note the recurrent presence of outliers at both tails of the various distributions totaling to more than 25% of the whole sample (N > 200). This indicates that the formation (or removal) of an IMHB can still have a significant impact on the biological activity independently from the IMHB features listed here. Importantly, owing to the symmetrical nature of the distributions, no significant enrichment in positive or negative outliers has been observed, implying that substantial increase and reduction in biological activity are equally probable. These results are consistent with earlier MMP-based analysis of substituent effects in medicinal chemistry ²¹⁹, highlighting the fact that bioactivity change outliers spans through different target classes, molecular topologies, transformations, and IMHB atomic pairs ²²⁰.

This analysis suggests that forming or disrupting an IMHB via minimal structural modifications will result in less than 0.5 log units change of biological activity more than half of the time. The general lack of directionality in the biological activity change and the occurrence of positive and negative outliers highlight the difficulty in generalizing the effect of IMHB on biological activity, as IMHB stem from specific electronic, conformational, and environmental constraints.

Despite the obvious limitations of inferring IMHBs based on molecular topologies and using biological data sources of likely different variability, the current data reinforce the notion that IMHB-mediated effects are highly context-dependent. Awareness of the underlying factors affecting IMHBs and the availability of case-specific structural information on the molecules of interest are essential prerequisites toward the fruitful exploitation of IMHBs in molecular design ²²⁰.

Studies of Silla el al. on the stereochemistry of α- and β-D-glucose stated that intramolecular hydrogen bonds are not responsible for the counter-clockwise arrangement of hydroxyl groups in D-glucose. The repulsion between the lone pairs of electrons in the endocyclic oxygen with oxygen bonded to the anomeric carbon dictates the preferred conformations of this compound, since this interaction induces the orientation of the remaining hydroxyl groups. These findings are instructive because they suggest that derivatization of D-glucose through replacement of hydroxyl by OR groups, which should not exhibit H-bond, can lead to counterclockwise arrangement, such as the prevalent rotamers of D-glucose ²²¹.

Several publications have reported beneficial effects of intramolecular hydrogen bonding on ligand-receptor binding and rationalized this with conformational restriction in which the small molecule substituents are favorably aligned with the protein pockets ^{222, 223, 224}. In these cases, the intramolecular hydrogen bond is crucial and removal of either hydrogen bond acceptor or donor results in a drastic loss of binding affinity. Rational design of internal hydrogen bonds for conformational preorganization was pursued in scaffold replacements for a diverse set of targets ^{225, 226, 227, 228, 229}. In addition to its effects on receptor binding, an improved property profile has been associated with molecules containing intramolecular hydrogen bonds. Increased brain penetration and pharmacological activity was observed for NK1a receptor antagonists and attributed to the higher apparent lipophilicity due to intramolecular hydrogen bonding ²³⁰. A team at Takeda has reported improved oral absorption and an excellent pharmacokinetic profile for luteinizing hormone-releasing hormone receptor antagonists when an intramolecular hydrogen bond was established ²³¹. Studies on cyclic peptides support the notion that the ability to form internal hydrogen bonds is critical for passive membrane permeability.

Geometric parameters of all small-ring intramolecular hydrogen bonds lie outside the typically observed ranges for unconstrained systems. Kuhn et al. ⁶⁷ illustrated this for carbonyl and nitrogen acceptors. A linear correlation between bond length and angle was observed in all cases. Intermolecular hydrogen bonds are close to linear, preferentially with angles greater than 150º. Only seven-membered ring intramolecular hydrogen bonds come close to these values. Six-membered rings have angles between 130º and 140º but the same distances as intermolecular hydrogen bonds. Five-membered rings have longer distances and smaller angles just within the cutoffs applied here (decreasing the lower cutoff to 75º changes the statistics only marginally by∼1%). On these grounds it may be argued that five-membered rings should not be regarded as classical hydrogen bonds but, more general, as favorable electrostatic interactions. We will still cover them here, since the conformation-property relationships for larger rings described in Kuhn et al. article are applicable to such systems as well.

Modulation of the electron-transfer capability is very important for the biological activity of quinones and hydroquinones. Among the interactions that play a central role in this issue, the formation of inter- or intramolecular hydrogen bonds in these molecules plays a key role ^{232, 233, 234}. A recent electrochemical study about quinones possessing intramolecular hydrogen bonds (IHBs) shows that this interaction stabilizes the anion radical structure, leading to a shift in reduction potentials toward less negative values when compared with quinones without IHBs ²³⁵. IHBs have shown appreciable effects on the antioxidant properties of hydroquinones and related phenols ^{236, 237}.

Martinez-Cifuentes et al. ²³⁸ studied intramolecular hydrogen bonds using MEP (molecular electrostatic potential) and NBO (Natural bond orbital) calculations, stating that the latter provide a better quantitative description of the strength of IHBs in o-carbonyl-hydroquinones, being more suitable to understand and predict the characteristics of this interaction. These results not only might be of interest to gain insight into intramolecular hydrogen bonds but also can help to rationalize the design of new hydroquinones with biological activity.

Studying the impact of stereospecific intramolecular hydrogen bonds on cell permeability and physicochemical properties, Over et al. ²³⁹ on series of stereoisomeric T. cruzi inhibitors, found an unexpectedly revealed striking influence of stereochemistry on solubility, lipophilicity, and cell permeability. This effect on compound properties was traced to the relative stereochemistry at two adjacent stereocenters where an intramolecular hydrogen bond was favored that reduced the basicity of the tertiary amine of the inhibitors while simultaneously shielding polarity from the surrounding environment. As a consequence, they were more lipophilic, less soluble, and had higher cell permeabilities than the stereoisomers for which an intramolecular hydrogen bond was less favorable. In comparison, limited profiling of other. Another unexpected finding was that the intramolecular hydrogen bond that influences the properties of this stereoisomers came via formation of an eight-membered ring. In contrast, the majority of intramolecular hydrogen bonds found in a recent exhaustive analysis of crystal structure databases involve formation of five- or six-membered rings ⁶⁷. Their observations emphasize the importance of preparing and screening pure stereoisomers in chemical probe or drug discovery projects, since their physicochemical as well as pharmacokinetic and pharmacodynamic properties may be significantly different.

These results obtained for the stereoisomeric T. cruzi inhibitors point to the opportunity that intramolecular hydrogen bonding can be used to “hide” hydrogen bond donors and adjust pKa in design of druglike compounds with properties at or beyond the RO5 (Lipinski’s rule of five). They also suggest that not only thermodynamically favored five- and six-membered rings ⁶⁷ but also intramolecular hydrogen bonding leading to formation of larger rings may be used in optimization of compound properties. It can be assumed that opportunities to adjust physicochemical properties, and subsequently cell permeability and oral bioavailability, will be of increasing importance as compound properties deviate further and further beyond the RO5. This is supported by the fact that intramolecular hydrogen bonding confers cell permeability and oral bioavailability to cyclosporine A, a cyclic undecapeptide ^{207, 240}, with properties far outside of the RO5. As revealed in this work and suggested by Alex et al. ²⁰⁷, conformational calculations may be used for prediction of intramolecular hydrogen bonding, and they may therefore have value as prospective tools for design and optimization of bRO5 compounds. Computed Δ log P_oct−tol values, i.e., differences in log P determined for partitioning between water and octanol or toluene, respectively, have also been suggested for prediction of intramolecular hydrogen bonding ²⁴¹. In addition, a recent analysis ⁶⁷ of crystal structure databases gave a list of intramolecular hydrogen bonding motifs for five- to seven-membered hydrogen-bonded rings that can also be deployed in compound design or optimization.

Modulation of challenging targets with extended binding sites requires compound classes that reach into chemical property space near the limit of what is acceptable for cell permeability and oral bioavailability ^{242, 243}, i.e., into bRO5 (beyond RO5) space ^{207, 244}. Macrocycles constitute one example of compounds that predominantly reside in bRO5 chemical space and that also have demonstrated success in modulation of challenging targets ²⁴⁵. A recent comprehensive investigation of macrocyclic drugs and clinical candidates revealed that a significant number that are orally bioavailable had molecular weights, lipophilicities, and polar surface areas that were higher than for traditional oral small molecule drugs ²⁴⁶. However, it was discovered that macrocycles, just as small molecule drugs, may have no more than five hydrogen bond donors to allow for oral administration ²⁴⁶. This observation, in combination with the findings reported herein, further emphasizes that masking of hydrogen bond donors by logical incorporation of intramolecular hydrogen bonds may be of particular value in efforts to improve cell permeability and oral bioavailability of compounds at the border of, or beyond, RO5 chemical space.

Previous data showing that several chloroquine analogs containing an intramolecular hydrogen bonding motif were potent against multidrug-resistant P. falciparum, led Madrida et al.to the exploration of the importance of this motif ²⁴⁷. So, 116 compounds containing four different alkyl linkers and various aromatic substitutions with hydrogen bond accepting capability were synthesized, as simple modifications, significantly altering the pKa and sterics of the basic side chain in chloroquine analogs. The obtained results may prove to be part of a strategy for overcoming the problem of worldwide resistance to affordable antimalarial drugs.

An attractive alternative to the permanent removal of HBD, to maintain target activity, is the introduction of a complementary HBA that is able to assume a conformation where the two groups can adopt a suitable geometry and distance to form a temporary ring or conformation where an intramolecular H-bond is formed. This temporary ring is different than the intermolecular conformations with solvent resulting in a thermodynamic equilibrium between the closed (intramolecular H-bond) and open (intermolecular H-bond) conformation. Intramolecular H-bonds, though often weaker than their intermolecular counterparts, have significant influence on properties such as charge distribution within molecules, the relative stability of conformers, and reactivity ²⁴⁸. The closed form, masking polarity from the environment and effectively removing one donor and one acceptor atom from the molecule, may then be less recognized as a P-gp substrate. Additionally, the closed form tends to be more lipophilic and might display faster membrane permeability ^{64, 207}. Indeed, Wright and Painter used uptake into red blood cells as a surrogate for membrane permeability to show that introduction of an intramolecular H-bond increased permeability of positional isomers of hydroxybenzoic acid ²⁴⁹. Despite potential ramifications for potency, removal of H-bonding groups is a much more effective and efficient strategy than adding more hydrocarbon and, in lieu of permanently removing H-bonds to maintain target activity, introduction of a intramolecular H-bond to mask groups that are required for potency may be the preferred strategy to achieve improved target engagement.

As weak intermolecular interactions play an important role in stabilizing a ligand energetically at the interface of a protein structure, the importance of hydrogen bonds in the binding affinity of a target-drug has been described extensively ²⁵⁰. But how hydrogen bonds optimize the hydrophobic interactions at the protein-ligand interface that increases the binding affinity of complex molecules, has not been properly analyzed ²⁵¹. For example, it has been done analyzing the hydrogen bonds and observed weak C-H…O interaction with an H…O distance of .1.85 Å. This kind of weak hydrogen bond can be broken and exchanged for another kind of bond, depending upon the chemical environment at the target, ligand and target-ligand interface.

To prevail in water environments, soluble proteins must protect their backbone hydrogen bonds from the disruptive effect of water attack by clustering nonpolar residues around them ^{252, 253, 254, 255, 256, 257, 258, 259, 260}. This exclusion of surrounding water, or wrapping effect, also enhances the electrostatic contribution by modulating the local dielectric (descreening the partial charges) and thus stabilizes the HB. Thus, underwrapped interactions, called dehydrons, represent thermodynamically unstable, vulnerable sites ^{252, 260} and the level of underwrapping has been shown to correlate with the degree of structural disorder ²⁵⁷ (in fact, the most extreme case is represented by the prion proteins which, being the worst wrapped proteins in the Protein Data Bank, are so thermodynamically unstable in monomeric form in solution that they aggregate to form amyloids ²⁵²; also among the worst wrapped proteins we can find the toxin proteins, whose dehydrons-rich structure can only be sustained by the establishment of disulfide bonds ²⁵²). A main notion that derives from this picture is that dehydrons are adhesive ^{252, 253}, hence promoters of molecular associations because their inherent stability increases upon approach of additional nonpolar groups. Thus, dehydrons constitute key motifs that signal protein binding sites. Under this scenario, the integrity of the interface of a biomolecular complex (when a protein binds to another protein or to a small-molecule ligand) becomes extremely reliant on intermolecular cooperativity based on three-body correlations ^{252, 253, 254, 255, 256, 257, 258, 259, 260}: a third nonpolar body protects an electrostatic interaction pairing the other two and not all three bodies belong to the same molecule. Thus, wrapper (nonpolar) groups become relevant for binding interactions when an interfacial hydrogen bond relies on them in order to remain over a critical wrapping value essential on stability terms. Such a decomposition of the complex interface into a web of three-body cooperative interactions enabled us to successfully predict the hot spots reported by alanine-scanning experimental studies for a set of protein-protein complexes ²⁵⁹.

The acidity ‘logK_a’ term for a given functional group relates to the strength of H-bond formed by the solute donor when it interacts with lone pairs of acceptor groups in solvent molecules. The basicity ‘logK_b’ term of a group relates to the strength of Hbond formed by the lone pairs of an acceptor (the solute) when it interacts with donor solvents ²⁶¹. Using such information, it is now possible to predict the relative strengths of individual H-bonds, or dipole–dipole attractions, and thus to quantify the various intermolecular forces between a drug molecule and its surrounding solvent molecules and/or receptors.

Several authors have reported surrogate in silico methods to calculate Abraham’s H-bond strength parameters to enable use of these parameters as descriptors for SAR analysis. Japertas et al. ²⁶² reported a fragment-based method while a more rigorous method based on quantum mechanics calculations has been described by Gancia et al. ²⁶³ The latter method was based on the premise that Abraham’s H-bond scales could be related to parameters derivable from quantum-mechanical calculations such as atom self-polarizability and electrophilic superdelocalizability. These are, broadly speaking, parameters that include terms relating to the degree to which a given molecular orbital receives a contribution from an atomic orbital, and the energy of the orbital. The calculations were performed using semi-empirical quantum mechanics methods such as Austin Model 1 (AM1) ²⁶⁴ or Parametrization Method 3 (PM3) ²⁶⁵. In some of the examples described in this review, we have used logK_a and logK_b values (as HBD/HBA strengths) calculated using an internal implementation of the quantum mechanical approach described by Gancia et al. ²⁶³

5. Tautomerism and binding

Tautomerism by definition concerns all molecules which can readily interconvert into isomers by transfer of a chemical group. (We refer to the commonly drawn structure of substances as “compounds” and to rare tautomeric form(s) of the compound as “tautomer(s).” Here “commonly drawn” means the form in which a molecule was saved in a database.) Tautomerism is very complex and is related to several phenomena: different types of migrating groups [an electrofuge or a nucleofuge], cationotropic and anionotropic properties, valence tautomerism, zwitterionic tautomerism, tautomerism related to migration of neutral groups in molecules, migration of bonds or ring-chain tautomerism ^{266, 267, 268}. As a subtype of cationotropy, the proton migration (prototropy) is the most commonly known tautomeric phenomenon. It concerns the movement of H atom between discrete sites of the ‘‘same’’ molecule (in contrast to ionization or protonation where H atom leaves or comes from another molecule). Even if only the case of prototropic tautomers is considered, the problem of tautomery in chemical collections remains complex.

One of the problems related to tautomeric representation is the prediction of molecular descriptors; since the outcome of every calculation crucially depends on what tautomeric form the calculation was based on. The following example shows the problem of log P prediction if the ‘‘wrong’’ tautomer is used and it is ignored, which tautomer in octanol/water partition may be favored. The molecules 4-hydroxypyridine (48) and 4-pyridone (49) (Figure 11) have a calculated CLOGP=0.93 for 48, and CLOGP=-1.31 for 49! ²⁶⁹ [CLOGP values varied based on method, nevertheless, the difference between 48 and 49 remains huge]. However, in solution tautomer 4-pyridone (49) is the predominant molecule. The experimentally measured log P for 4-hydroxypyridine is=-1.3, which is very close to the value calculated for 49. This clearly illustrates that obtaining the appropriate form of the tautomer improves the quality of the property prediction.

The last example concerns the intramolecular H-bond due to the tautomeric interchange ^{270, 271}. López-Rodríguez et al. ²⁷¹ studied the prototropic equilibrium in the series of serotonin 5-HT4 receptor ligands. Two classes of benzimidazole derivatives were analyzed to gain insight into the bioactive conformation of these novel ligands. Their results from NMR and IR techniques and theoretical methods confirm the presence of important intramolecular hydrogen bonds between the benzimidazole ring and adjacent side chain groups. These hydrogen bonds are possible also due to the tautomeric hydrogen shift on benzimidazole ring. Thus, the molecular skeleton as well as the energy required for conformational changes has the effect on adopting a bioactive conformation for ligand–receptor interaction. Such structural studies [see also in Ref. ²⁷¹] provide important insights to guide the design and synthesis of new compounds with predetermined pharmacological activities.

The tautomeric equilibrium is influenced by a number of variables including concentration, temperature, pressure, solvent type, and pH. Tautomers differ in heat and energy of formation, proton affinities, dipole moments, and ionization potentials [see examples in Refs. ^{272, 273}]. There are numerous studies of tautomers in gas or aqueous solution, however, little is known about tautomerism of ligands in the binding site of proteins. Let us imagine a simple case of a transfer of one hydrogen atom on the same molecule-ligand. Tautomery can alter the skeleton of a given molecule, which in principle can be seen as a new distinct molecule with different complementarity to the target. Thus, Sanchez-Moreno et al. studied the differences in activity against T. cruzi Fe-SOD enzyme of the possible tautomers of several imidazole derivatives ^{274, 275, 276, 277, 278}.

Thus, several questions arise: does a molecule bind preferably in one distinct tautomer? Is the most stable tautomeric form in aqueous solution also the most stable form in the active site of a protein? What can be the binding contribution of a ligand in its excited tautomeric state in contrast to its ‘‘normal’’ tautomeric state, e.g., its low energy configuration? How to treat a compound, whose proton shift induces different stereoisomery? Compared to the large amount of data available on ligand-binding interactions, little is known about the binding modes of distinct tautomers of a ligand molecule.

The general view is that the binding environment within a protein is a very specific one. It is different from the environment of the aqueous solution or the vacuum. Apolar or polar, acidic or basic side-chains create local pHs, shift their pK values and subsequently influence the functional groups of the ligand. Presence of a ligand, metal cations, and water also influence the pH (and the pK) conditions in close proximity of amino acid side-chains and the process of catalysis ²⁷⁹. In such a context, ligands may be ionized or can achieve its excited tautomeric state. Tautomers have different molecular shapes and different hydrogen-bond donor and acceptor properties resulting in a significant impact on molecular recognition. As it was shown in this section, tautomers are able to bind to the active site of a protein. Nevertheless up to now, tautomers have been neglected, even omitted in automated molecular docking.

Recent advances in molecular docking occurred in various fields such as ligand–protein flexibility, scoring function or automatized processing ^{280, 281, 282, 283}. However, little progress has been achieved in the simulation of the immediate environment within the active site. This would demand predictions of terms such as solvation-desolvation, temperature of the system, microenvironment of the active site in proximity of the ligand and protein side-chains, pH, ionization, protonation/deprotonation, and tautomerism of ligands upon binding.

Including tautomerism in virtual screening procedures should improve the reliability of the screening due to the following reasons: (i) it enlarges the chemical space covered by the database and (ii) it takes into account that a compound can bind in its tautomeric state, which raises the chance of detecting a hit. Hence, this approach can be considered as a large improvement of the virtual screening procedure itself. Subsequent selection of the best-ranked compounds or tautomers may also reveal new characteristics of ligand binding. If a docked compound is stabilized in its tautomeric state in a given binding site environment, this can lead to the recognition of a putative ligand, which would not be detected using classical screening protocols.

Few investigations have addressed this topic. ProtoPlex, a program developed by Pearlman et al. generates all protomers of drug-sized compounds (tautomer and/or protonation state) and yields false negatives if the protomer in the screening library is not preferred by the receptor ²⁸⁴. Similarly, Sadowski et al. addressed the issue of tautomery together with protonated molecules using a tautomer and protonation preprocessor for virtual screening ²⁸⁵. Pospisil et al. have used the generator of tautomers AGENT to generate a database of energetically probable tautomers, which in parallel with the database of compounds can be docked to the target of interest ²⁸⁶.

6. Calculation of Hydrogen Bonds Characteristics

Because of the importance of hydrogen bonds to drug design, much work has been done in the past on the theoretical modeling of hydrogen bonds for QSAR studies. The approaches span from simple indicator methods ^{287, 288} to parametrization approaches using theoretically calculated properties, such as atomic charges ²⁸⁹, molecular electrostatic potential ²⁹⁰, LUMO and HOMO properties ^{289, 291}, atom polarizability ²⁹², and superdelocalizability ^{263, 293}, to model hydrogen bond strength. However, due to the character of being the property of a group of atoms and the susceptibility to local environments, hydrogen bonds could not be modeled accurately by a general semiempirical or rule-based method because there are many exceptions, such as steric factors, to be accommodated by a finite set of rules. To treat complicated systems such as hydrogen bonding, the ab initio quantum chemistry method is an ideal approach since all electronic and steric effects are fully taken into account in such a treatment. Over the past decade, much work has been published on the ab initio study of hydrogen bonds ^{70, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307}. However, few ab initio calculations have been carried out on hydrogen bond systems that are directly related to medicinal chemistry in drug design.

Approaches for estimating the strengths of H-bonds can be divided into two groups, quantum chemical calculations and empirical methods. The use of ab initio calculations allows the energy and geometric structures of small molecules to be described with an accuracy approaching those of experimental results. However, a similar treatment of large molecules is not yet practical. Reliable empirical H-bond parameters for drug-like molecules need to be derived from large databases of experimental measurements.

Raevsky’s program HYBOT is linked to a database of several thousand thermodynamic measurements of 1:1 H-bonding systems and can estimate H-bond donor and acceptor factors for the relevant sites in any compound based on its chemical structure ³⁰⁸. There is however a difference between H-bonding scales that refer to 1:1 complexation and scales that refer to the overall or summation H-bond acidity and basicity. The second scales deal with the situation in which a solute is surrounded by solvent molecules so that all acid (i.e. H-bond donor) groups and all basic (i.e. H-bond acceptor) groups in the solute are potentially involved in acid–base interactions at the same time. In the approach by Abraham ³⁰⁹ these differences are explicitly addressed and it has been concluded ³¹⁰ that functional group constants derived from overall H-bond parameters are the more useful in analyses of physicochemical and biochemical properties.

Property-based design ³¹¹ is a difficult problem that requires, on the one hand, an understanding of the properties of biological barriers and, on the other hand, an adequate and consistent description of the properties of drug molecules. Simple parameters (e.g. number of H-bond donor and acceptor sites) are highly useful but we believe that a more rigorous approach is desirable for the description of molecules and advocate the use of Abraham descriptors and the general solvation equation. In the following sections this approach is explained in more detail and some applications are examined.

Hao theoretically studied hydrogen-bonding strength for drug molecules using quantum chemical calculations carried out with the program package Jaguar from Schrodinger, Inc. ³¹² The standard energy difference method ^{70, 313} was used to calculate the energy of the hydrogen-bonding interaction between a donor and acceptor

(10)

where ΔE_HB is the energy of bonding interaction, R is the set coordinates that define the structure of the hydrogen-bonding complex, E_AD(R) is the total energy of the complex, and E_A and E_D are the individual energies of the donor and the acceptor, respectively.

One of the major technical issues encountered in his work was to find a proper procedure for geometry optimization. For simpler hydrogen bond systems, such as when one subunit is a water molecule, a single constraint on the hydrogen bond angle would be sufficient to keep the donor and the acceptor in a proper geometry through the geometry optimization process ³¹⁴. However, for calculations involving relatively large donors and acceptors in which only one or two internal variable constraints were applied, geometry optimization by energy minimization often results in structures that are substantially different from the desired hydrogen bond conformation. The reason for such an outcome is that, when two relatively large molecules are free to move, energy minimization always brings secondary interactions, in addition to hydrogen bonding, into the total energy. The secondary interactions not only perturb the hydrogen bond geometry but also cause difficulty in separating hydrogen-bonding energy from other nonbonded contributions.

In a realistic and complete inhibitor-receptor complex, the hydrogen bond interactions between the donor and the acceptor are not totally free to reach an absolute energy minimum but are subjected to numerous additional constraints at the binding site such as tethering from the protein structure. As a result, the hydrogen bond geometry in inhibitor-receptor complexes can be quite different from that produced from energy minimization of free subunits. We desired a procedure that computes hydrogen-bonding interaction at a specified geometry and can control the effects of secondary interaction contributions to hydrogen-bonding energy.

Hao demonstrated in this work that ab initio calculated hydrogen-binding energy has an excellent linear correlation with logarithm of experimental hydrogen-bonding constants. This provides a basis for the theoretical prediction of hydrogen-bonding strength for a series of acceptors (or donors) with respect to a given donor (or acceptor). The ab initio approach advances the level of sophistication of theoretical modeling of hydrogen bonds and may provide deeper insights to the mechanism of hydrogen-bonding interactions. The method also provides a tool for handling steric effects, conformational properties, and secondary interactions of specific systems in predicting hydrogen bond strengths. These effects were traditionally not treated in empirical and parametrization approaches for hydrogen bonds and so he stated that his method was helpful in dealing with hydrogen bond related SAR problems in drug design.

Another example of computational involvement in drug design is molecular docking. The aim of molecular docking is to optimize the geometry of both the receptor, the drug and their orientation so that the Gibbs free energy of the system is minimized ^{315, 316}. Docking assesses how well the compound fits into the active site of the receptor based on shape complementarity, conformation and the energy of interactions made with the receptor. The molecule with the best fit is then optimized to ensure that all other possible binding orientations have been considered as more than one conformation of the molecule may fit into the active site. The docked molecule is then simultaneously further analyzed by adding or removing functionality, the program can then score how this affects the docking score and flag any potential lead compounds for further investigation.

The receptor and ligand can adopt many conformations. In order to conserve computational time most docking searches consider the protein as structurally rigid and the ligand has full conformational freedom. This means that some binding modes where there is a different protein conformation maybe missed. Docking algorithms such as GOLD (Genetic Optimization for Ligand Docking) use a genetic algorithm to screen more binding modes by allowing flexibility of the protein active site. This effectively allows the structure of the protein active site to simultaneously evolve as the conformation of the ligand is changed. GOLD also allows the user to consider water-mediated binding between the drug and the receptor ³¹⁷.

The investigations of structural conformers, molecular interactions and vibrational characterization of pharmaceutical drug are helpful to understand their behavior. In the work of Karthick et al. ³¹⁸, the 2D potential energy surface (PES) scan was performed on two dihedral angles to find the stable conformers of busulfan. The strong and weak hydrogen bonds between the functional groups of busulfan were then analyzed using quantum topological atoms in molecules analysis, giving AIM analysis the topological picture of hydrogen bonds and the ring formation. A couple of cage critical points were found in some conformers that may cause their predominant stability. With the help of RDG plot and isodensity surface analysis, the interactions taking place within the molecules was classified.

The most common approach for the investigation of the energetics of intramolecular H bonds is usually based on analysis performed by ab initio methods ^{319, 320, 321, 322}. Generally, the strength of the hydrogen bond is calculated as the difference between the enthalpies of the structures with and without this bond.

The isodesmic reactions appear to provide more reasonable results than the cis−trans approach when used to evaluate the intramolecular interaction energy ^{323, 324, 325}. In isodesmic reactions ¹⁸⁷, the number and types of bonds are conserved on the reactant and product sides of the reaction. The isodesmic method was used for the estimation of intramolecular HB energies based on the assumption that the total molecular energy can be partitioned into energies of chemically recognizable fragments ^{325, 326, 327, 328, 329}. However, the isodesmic reaction approach does not provide the true H bond energy but includes the effect of the energy due to the formation of the ring structure ³³⁰. This method is advocated for systems with one HB but is not recommended for the estimation of the single intramolecular H bond energy in polyhydroxy systems ³²⁷.

Hydrogen bond is an important type of interaction between drug molecules and their receptors, and so it has been presented a computational method for accurately predicting the hydrogen-bonding strength for different acceptors with respect to a given donor or vice versa. The method is based on quantum chemistry DFT calculation of the interaction energy between hydrogen bond donors and acceptors. An excellent linear correlation is observed between the calculated hydrogen-bonding energies and the experimentally measured hydrogen-bonding constants log K_β on a variety of known hydrogen bond acceptors and donors. These results not only indicate the predictive power of this method but also shed light on factors that determine the magnitude of experimentally measured hydrogen-bonding constants for different acceptors with respect to a given donor, suggesting a primarily enthalpic contribution from hydrogen-bonding energy. The method can be used for evaluating the effects of steric interference and inhibitor binding geometry on hydrogen-bonding strength in drug design ¹⁸⁴.

This novel method proposed for the estimation of intramolecular hydrogen bond energy is based on fragmentation strategies. This new simple computational method uses a density functional theory (for preliminary optimization) and an ab initio procedure during final optimizations and single point energy calculations. The correctness of the modification of the MTA (molecular tailoring approach) method was verified using a number of representative intermolecular interactions. This fragmentation method was tested on over 100 hydroxycarbonyl compounds wherein only three single point energies of tailored fragments are required for the estimation of one intramolecular hydrogen bond energy. The estimated intramolecular hydrogen bond energies range from −1.4 to −13.7 kcal/mol and show a qualitative rank ordering with the O−H···O corresponding lengths, distance d_(O···O), stretching frequencies v_OH, NMR chemical shift δ_H, and electron density topological parameters ρ_BCP and ∇²ρ_BCP at a (3,−1) value. Moreover, it appears that for this group of compounds the energy of intramolecular hydrogen bonding is not bound with the electron density in the hydrogen-bonded ring critical point.

The values of the hydrogen bond energy calculated by this MTA-modified method are not always in agreement with those reported in the literature, although they are generally consistent with the parameters typically used to describe the problem. The model is able to convincingly interpret the fine changing of structural fragments as a cause of the weakening/strengthening of the intramolecular hydrogen bond by electrondonating/withdrawing substituents or conjugated extra hydrogen bonds accepted by the carbonyl. The hydrogen bond strength for all of the compounds depends on the spatial arrangement of the bonds, steric accessibility of the donor−acceptor environment, and cooperativity/anticooperativity with other HB bonds. Importantly, even the subtle structural and stereoelectronic effects are well reflected by the hydrogen bond relative energies. Moreover, Rusinska-Roszak et al. ¹⁸⁴ show that when the molecules become more complex or the number of intramolecular hydrogen bonds increases (including the cooperative H bond phenomenon) the presented exhaustive stereoelectronic effects interpretation may still be accessible.

The H bond as an intrinsic feature of the ground-state structure of many molecules may affect their shape and properties. The calculated values of E_HB can be a good and intelligible explanation for the number of interactions and the structural and conformational H-bonding phenomenon in a given compound. The Rusinska-Roszak et al. method, being quite general, can be applied also to the resonance-assisted hydrogen bonding, to the aromatic systems and, generally, to intramolecular hydrogen bonded systems, including N−H···O and O−H···N in amides and peptides. The five-membered structures of the O−H···O=C intramolecular hydrogen bond are essentially dependent on the ring molecular strains and special steric interactions.

In Nocker et al. work ³³¹ hydrogen bonds are calculated by development of a more general method avoiding lots of restraints basing on the supramolecular approach and use of Density Functional Theory (DFT), a well-established method when dealing with this type of interactions ^{332, 333, 334, 335, 336, 337, 338}, that in comparison with force-field derived hydrogen bond energies points out an underestimation of those in all studied force fields ³³⁹. DFT expresses the electronic energy of a system in terms of its density ³⁴⁰ instead of using the many-electronic wave function and depends on the knowledge of the exchange correlation energy functional for which different approximations can be used ³⁴¹. Interaction energies between structurally different acceptors and a common hydrogen bond donor as well as between different donors and a chosen hydrogen bond acceptor were examined by the use of the hybrid functional B3LYP and basis set aug-cc-pVDZ. Further calculations for investigating possible outliers were performed with the smaller 6-31+g* basis set and B3LYP as well as MPW1PW91 procedure. For comparison multilevel G3MP2B3 calculations with a high accuracy were performed, and experimental data were obtained from Abraham et al. ²⁵.

It is worth to note that obtained thermodynamic constants of complex formation are free energies, while quantum-mechanical studies lead to enthalpies. A linear correlation of the two quantities implies that entropy changes are constant on hydrogen bond formation or are at least linearly depending on its strengths. Influence of electron-withdrawing and -accepting substituents on the aromatic systems phenol as H-bond donor and pyridine as H-bond acceptor was elucidated by Nocker et al. ³³¹ using correlation with σ-Hammett constants.

Theoretical determination of hydrogen bond energies in Nocker et al. work ³³¹ is achieved by equations 10 (uncorrected) and 11

(11)

where ΔE_HB is the interaction energy of bonding, E_AD is the energy of the geometry optimized complex, and E_D as well as E_A are the energy of the individual geometry optimized hydrogen bond donor and acceptor.

In quantum chemical calculations correction of the Zero Point Energy (ZPE) and Basis Set Superposition Error (BSSE) is thought to be needed ³⁴². The ZPE summarizes the electronic ground state energy for each nucleus’ vibration considering them as harmonic oscillator. The BSSE is a systematic error when dealing with complexes formed by fragments ³⁴². It can be estimated with the so-called Counterpoise correction by Boys et al. ^{313, 343} that eliminates the error that arises from calculating complex and fragments with different basis sets.

Computational methods such as linear free energy relationships (LFERs) offer a useful high-throughput solution to quickly evaluate drug developability, e.g. membrane permeability, organic solvent/water partition coefficients, and solubility. LFERs typically assume the contribution of structural components/functional groups to the overall properties of a given molecule to be constant and independent.

Significant nonadditivity was observed in peptides in that the contribution of the peptide backbone amide to the apparent transfer free energy from water into the bilayer barrier domain is considerably smaller than that of a “well-isolated” amide and greatly affected by adjacent polar substituents on the C-termini. In order to explain the phenomenon of nonadditivity, the formation of intramolecular HBs and inductive effects of neighboring polar groups on backbone amide, were investigated using FTIR and MD simulations ³⁴⁴. Both spectroscopic and computational results provided supportive evidence for the hypothesis that the formation of intramolecular HBs in peptides is the main reason for the observed nonadditivity of Δ(ΔG°)-CONH-. The MD simulation results of Cao ³⁴⁴ showed that the inductive effect of neighboring groups is not as important as the effect of intramolecular HBs.

Because it is the properties conferred by the functional groups that are most important to the biochemical activity of a given compound, many CADD applications treat functional groups containing different atoms but conferring the same properties as similar or even identical. For example, the capacity for hydrogen bonding can heavily influence a molecule’s properties. These interactions frequently occur between a hydrogen atom and an electron donor such as oxygen or nitrogen. Hydrogen bonding interactions influence the electron distribution of neighboring atoms and the site’s reactivity, making it an important functional property for therapeutic design. Commonly, hydrogen bonding groups are separated as hydrogen bond donors with strong electron-withdrawing substituents (OH, NH, SH, and CH) and hydrogen bond acceptor groups (PO, SO, CO, N, O, and S) ^{345, 346}.

Consider the following hydrogen bond equilibrium between the hydrogen bond donor 4-fluorophenol (blue) and the hydrogen bond acceptor acetone (red) (Figure 12).

The equilibrium constant K_HB for the hydrogen bond acceptor can be measured and converted into pK_HB (equations 12-13) and also for the hydrogen bond donor (equations 14-15). Increasing the polarity of either or both the donor and acceptor will increase the pK_HB value.

(13)

(14)

(15)

In Figure 13 the two chlorines are identical, which means they can both interact with the hydrogen bond donor. When measuring the pK_HB values for compounds containing more than one identical acceptor or donor functional group statistical corrections have to be made as an interacting molecule will bind to both chlorines giving a larger K_HB, normally a log_n correction (where n is number of identical sites) is applied.

The first pK_HB scale was devised by Taft et al. ³⁴⁷. Their studies involved ¹⁹F NMR (Nuclear Magnetic Resonance) titrations. This type of technique observes the fluorine nucleus being deshielded upon binding giving rise to a change in chemical shift (Δδ) between the bound and unbound 4-fluorophenol. The Δδ is measured for different concentrations of acceptor until saturation occurs. This data can then be fitted to obtain K_HB. These studies required the use of an internal/external reference standard which was 4-fluoroanisole as this closely resembled 4-fluorophenol without the hydroxyl proton to make a hydrogen bond. This allows calibration for the hydrogen bond equilibria studied and proved to be a suitable internal reference even for weak bases where higher concentrations were required to achieve saturation.

One method for obtaining pK_HB values is via IR (Infrared) studies ^{348, 349, 350, 351, 352}. This is done by measuring the difference in IR stretch (ΔνOH), between the free and bound donor. The greater the ΔνOH, the higher the pK_HB. In these studies it was also shown that there is a functional group dependent correlation between ΔνOH and pK_HB) (equation 16).

(16)

A = gradient co-efficient, B = offset on pK_HB axis.

This is particularly advantageous as additional pK_HB values can be estimated for a donor that has poor solubility. pK_HB values obtained in this way are secondary values and the accuracy depends on the quality of the correlation between the ΔνOH and pK_HB for that particular functional group.

Many of the studies in the literature report cases where electronic effects predominate over steric effect of a molecule. Resonance effects were observed in a study on amides ³⁵⁰. Tertiary amides have a higher pK_HB than secondary amides due to more electron releasing alkyl substituents. A back donation effect is also observed when attaching phenyl rings to the nitrogen. Increasing phenyl substitution leads to an increased conjugation of the nitrogen lone pair (Figure 14).

A study on para and meta substituted acetophenones showed that in para substituted acetophenones the main contributor to pK_HB was π-conjugation ³⁵², this effect being slightly smaller in meta substituted systems. In both cases the substituents can also donate electron density through the σ-bonding framework inductively, but for para substituents resonance effects predominate. But in meta substitution these inductive effects are more profound. Extending the conjugation of a system leads to an increase or decrease in pK_HB depending on whether EDGs (Electron Donating Groups) or EWGs (Electron Withdrawing Groups) are attached.

Intramolecular hydrogen bonding can make a significant difference in the hydrogen bond donating and accepting strength of a molecule ³⁵³ In phenyl ester 59 (Figure 15), the two carbonyl lone pairs are available for making a hydrogen bond. Upon substituting the ortho position with a hydroxyl group, an intramolecular hydrogen bond in phenyl ester 60 increases the electrostatic potential around the carbonyl oxygen hence reducing the hydrogen bond acceptor ability of the carbonyl. Replacing this hydroxyl in phenyl ester 61 removes this intramolecular hydrogen bond, while the methoxy group is also electron donating via π-conjugation, so the basicity of the carbonyl is increased.

For hydrogen bond donors, the same steric argument as for acceptors applies. Larger groups surrounding the hydrogen will weaken the hydrogen bond.

To better define the nature of H-bonds in biological systems, two key aspects need to be considered. First, the protein-ligand H-bonding process will invariably compete with bulk water. This process is represented in Eq. 1, which shows a reversible competing H-bonding process between a ligand H-bond acceptor (L; note that L stands for a ligand atom, not the whole ligand) and a protein receptor H-bond donor (PH), where the square bracket indicates that the water molecule’s orientation and position is constrained by the protein

(17)

(18)

(19)

The free energy change (ΔG) of this process is usually considered as the Gibbs free energy change from PH···[OH₂] to PH···OH₂ (Eq. 18). The difference between Eqs. 17 and 18 describes an H-bond competing process associated with small molecules in aqueous solution (Eq. 19). Because this relationship is poorly understood and often ignored, a common misperception—that is, generating stronger protein-ligand H-bonds leads to a net gain in binding affinity—exists. The ΔG in Eq. 18 is not dependent on the strength of protein-ligand interactions, whereas the ΔG in Eq. 19 is associated with protein-ligand H-bonds. Thus, the ΔG in Eq. 19 provides useful quantitative information in deciphering how protein-ligand H-bonds may modulate ligand binding affinity. To address the first issue of competing H-bonds in bulk water, Chen et al. ³⁵⁴ proposed a new H-bond pairing principle to evaluate the ΔG in Eq. 19, and demonstrate that the nature of these H-bonds depends on the pairing of the donors and acceptors.

Second, H-bonding in biological systems is highly complicated. Some key determinants, such as solvent entropy changes during the H-bonding process, are difficult to measure accurately using either experimental or theoretical methods. This limitation is a significant reason why the contribution of H-bonds to biological function remains poorly defined. Furthermore, the net free energy contribution of an H-bonding process represents the sum of several parts, with individual values being much larger than the net contribution in some cases. Even if each component can be measured with small relative error, the net contribution cannot be obtained with accuracy. To address this second issue, we developed a novel parameter derived from experimental partition coefficients to calculate the contribution of specific H-bonds to ligand binding affinity. Because this parameter includes the factors that influence the free energy contribution of H-bonds, notably electrostatic interactions, desolvation, entropy change of solvent, and van der Waals interactions, this makes the calculation simple and accurate because summarizing the individual parts that are hard to quantify accurately is not necessary. By applying both the H-bond pairing principle and the novel parameter, we examined the mechanism and the extent to which protein-ligand H-bonds modulate ligand binding affinity.

The H-bond competing process can be defined by the following general equation, where two acceptors (A1 and A2) and donors (D1 and D2) form mixed pairings

(20)

When A₂ and H-D₂ have stronger H-bonding capabilities than A₁ and H-D₁, respectively, Eq. 20 favors (both in enthalpy and in free energy) the pairing A₂···H-D₂. Chen et al. ³⁵⁴ estimated the H-bonding capability of an atom using the free energy required to transfer the atom from water to hexadecane. They then used a modification of the method of Kenny et al. ³⁵⁵ to evaluate the H-bonding capability of the respective atoms in the protein-ligand complex by calculating the difference (ΔlogP₁₆) between water/hexadecane partition coefficients (logP₁₆) and logP₁₆ of saturated hydrocarbon molecules with the same molecular surface area.

It is known that ΔlogP between hexadecane (or alkane)/water and 1-octanol/water provides a measure of the H-bond potential of a molecule ^{356, 357, 358}, and the calculated standard errors for the H-bonding capabilities calculated in this way appear acceptable. The H-bonding capabilities of some atoms are highly accurate [low root mean square error (RMSE) values] because values can be obtained directly from one-step calculation (Eq.18) and the RMSE value for the basic equation (Eq. 17) is lower. For atoms where the experimental hexadecane/water partition coefficients for their relevant compounds are lacking, RMSE values increase (>1) because the calculation steps increase. Systematic measurement of experimental hexadecane/water partition coefficients for the compounds containing the atom type and/or the relevant compounds will improve future calculations of H-bonding capability.

Experimental support for the strong-strong or weak-weak (s-s/w-w) H-bond pairing principle is provided by the observation that the H-bond competing process between two halide ions and two hydrogen halides in a gaseous system (Eq. 21) favors the s-s/w-w pairing in enthalpy.

(21)

The enthalpy change (ΔH) in Eq. 21, calculated from the experimental data for the H-bond energies, is −19.9 kcal/mol ³⁵⁹, indicating that the equation favors s-s/w-w H-bond pairing in enthalpy because HCl is a stronger H^-bond donor than HF and F⁻ is a stronger H-bond acceptor than Cl⁻. This phenomenon is universal because all H-bond competing processes usually favor s-s/w-w pairing in enthalpy. Although entropy-enthalpy compensation exists, the favorable enthalpic contributions of H-bond competing processes are only partially canceled by unfavorable entropic contribution (TΔS). Thus, H-bond competing processes favor s-s/w-w pairing in free energy.

The ΔG for the reversible competing protein-ligand H-bonding process shown in Eq. 17 has two components: (i) the ΔG associated with the release of a well-ordered water molecule into the bulk solvent (Eq. 18), which does not depend on protein-ligand interactions, and (ii) the ΔG associated with protein-ligand H-bonds (Eq. 19). The ΔG in Eq. 19 cannot be obtained from experimental data. However, because the H-bond competing process between the same H-bonding protein atom and different ligand atoms obeys the H-bond pairing principle, ΔG can be calculated by comparing the experimental binding affinities of the two ligands. The ΔG for the H-bond competing process of two ligand atoms with the same protein atom(s) can be expressed as shown in Eq. 22 ³⁵⁴

(22)

where k is a constant and is equal to 1/H_WH; H_PH, H_WH, H_A, and H_B are the H-bonding capabilities of the protein atom(s), the H-bond donor or H-bond donor of water, and the atoms of ligands A and B, respectively. The effect of H-bond geometry on ΔG is factored into the H-bonding capability of the protein atom(s).

Equation 22 is a mathematical expression of the H-bond pairing principle. Because this derivation is complex, some approximations are used to develop the models. For example, to derive Eq. 22, it is assumed that single H-bonds of similar distance make up the pairings. However, the calculated ΔG for the H-bond pairing process—in which two H-bonding acceptors compete and bind to the same nonpolar site on a protein receptor—is the same as the ΔG obtained from the experimental water/hexadecane partition coefficients for any H-bond pairings. These experimental findings validate the model as fitting beyond its approximations. Chen et al. ³⁵⁴ further validate Eq. 22 by showing that strong H-bonds in different pairing models have opposing effects on experimental binding free energy. Their work shows that strong s-s pairing H-bonds between 63 and its protein receptor enhance binding affinity, whereas the strong-weak (s-w) pairing between 65 and its protein receptor is not as favorable as the w-w pairing provided in the form of a polar-apolar interaction. Moreover, the reported binding affinities of two structurally similar scytalone dehydratase inhibitors 62 and 63 ³⁶⁰ indicate that substitution of an apolar H atom (H-bonding capability, 0) for a cyano group (H-bonding capability, 16.0) enhances receptor antagonism by ~30,000-fold. The binding free energy difference (ΔΔG) between 62 and scytalone dehydratase inhibitors mainly results from (i) H-bond interactions with Tyr30 and Tyr50 (ΔG_HB) and (ii) the relative flexibility of the OH groups in Tyr30 and Tyr50 (ΔGflex) because the OH groups interacting with 63 are less flexible. The first term, ΔG_HB, is the ΔG of the H-bond competing process, which is −33.2 ± 3.2 kJ/mol because the H_PH, H_A, and H_B of the process are 21.6 ± 1.5, 0, and 16.0 ± 0.5, respectively.

The free energy required for fixing two rotatable bonds (ΔGflex) is ~7.4 ± 1.8 kJ/mol because the predicted free energy cost for rotor restrictions is close to 3.7 ± 0.9 kJ/mol per rotor ^{361, 362}. Thus, the calculated ΔΔG is −25.8 ± 3.7 kJ/mol, which is in close agreement to the experimental ΔΔG (−25.6 kJ/mol at 298 K). On this basis, Chen et al. ³⁵⁴ conclude that the s-s pairing H-bonds between the cyano group of 63 and the receptor tyrosine hydroxyls can markedly increase the binding affinity. Further evidence that shows that the H-bond interactions between the inhibitor 63 CN group and the receptor tyrosine hydroxyls are strong and favorable to binding affinity is based on their geometry and large effects on binding affinity.

By contrast, the reported binding affinities of two heterocyclic aromatic sulfonamide inhibitors of carbonic anhydrase (64 and 65) indicate that the strong H-bond between 65 and Thr200 is not as favorable as the weak (polar-apolar) interaction between 64 and Thr200. The binding affinity of 65 is ~30-fold lower than 64. Two factors contribute to the difference in binding affinity: (i) differential H-bond interactions with Thr200 (ΔG_HB), which is equal to the ΔG for the H-bond pairing, and (ii) the free energy difference in transferring the hydrogen atoms from water to protein (ΔG_sol). Because protein has a higher dielectric constant than hexadecane, ΔG_sol is less than the free energy difference in transferring the hydrogen atoms from water to hexadecane (which is 8.4 kJ/mol based on the H-bonding capabilities of the hydrogen atoms). Thus, transferring the hydrogen atom of 65 reduces the activity <30-fold (=108.4/5.71) compared with the hydrogen atom of 64. Because the difference in activity between 65 and 64 is ~30-fold, the strong H-bond between 65 and Thr200 (s-w pairing) is not as favorable as the weak H-bond between 64 and Thr200 (w-w pairing).

Concordance between experimental and calculated results is supportive of Eq. 22, and from this, the free energy contribution of a protein ligand H-bond (ΔG in Eq. 19) can be calculated from Eq. 23

(23)

where k is a constant and H_W, H_L, and H_P are the H-bonding capabilities of the H-bond donor/acceptor of water, ligand atom, and protein atom, respectively. The relationship between ΔG_HB and H_L, H_P offers insight into how specific protein-ligand H-bonds contribute to binding affinity in the following ways:

(i) H-bonds with s-s pairings increase ligand binding affinity. An example is provided by the strong binding affinity of the scytalone dehydratase inhibitor 63 ³⁶⁰. The s-s H-bond pairings between the cyano group (H_L = 16.0 > H_W) and two strong donors from Tyr30 and Tyr50 (H_L = 21.6 > HW) increase the binding affinity by ~300-fold. Another noteworthy analogy exists in nature, where the high-affinity binding of biotin to streptavidin (Ka = 2.5 × 1013 M−1) ³⁶³ represents a prototypic example of how s-s H-bond pairings facilitate exceptionally strong interactions for a molecule with such few heavy atoms. All H-bond–forming atoms from the ureido ring of biotin form s-s receptor-ligand H-bond pairings. The H-bonds for the ureido oxygen contribute significantly to the binding affinity because of the extreme H-bonding capabilities of both donors and acceptors. For the H-bond acceptor, the ureido oxygen has two lone pairs of electrons with an H-bonding capability of 14.3 per lone pair. For the H-bond donors, the sum of the H-bond capability of the three H atoms in Asn23, Ser27, and Tyr43 is 25.0, much stronger than the sum of the H-bonding capability of two water H-bond donors (14.04). Critical importance of this s-s H-bond pairing in reducing competitive interference with water is demonstrated by the biotin analog 2-iminobiotin, where binding affinity to streptavidin is reduced by more than 3 million–fold ³⁶⁴. Although there is some unfavorable positive-positive repulsion between 2-iminobiotin and the side chains of streptavidin, these interactions are minor because these side chains are rotatable. Thus, the s-s H-bond pairings of the ureido oxygen atom are important for the high affinity of biotin to streptavidin.

(ii) Ligand binding affinity is relatively unaffected by the strength of H-bonds when the H-bonding capability of the ligand or protein atom is close to that of water (ΔG_HB ≈ 0). The binding of the heterocyclic aromatic sulfonamide inhibitors (66, 67, and 68 in Figure 17) to carbonic anhydrase ³⁶⁵ represents an example where significant modification of H-bond strength fails to improve the binding affinity. The ligand acceptors form H-bonds with OH and NH of Thr200, but the latter H-bond is weak because of large donor-acceptor distances. Because the H-bonding capability of the protein is close to that of water, ligands 66, 67, and 68 have similar binding affinities, although the H-bond acceptors have markedly different H-bonding capabilities.

(iii) s-w H-bond pairings are less favorable than w-w pairings, although they might be stronger. A typical example in Figure 18 indicates that the nature of the H-bond pairing is more important than H-bond strength. This phenomenon is further illustrated by the favorable interactions between aromatic rings and polarized CH groups ^{366, 367}. Aromatic rings are weak H-bond acceptors ^{366, 368}, whereas CH groups are weak H-bond donors ^{368, 369}. The favorable hydrophobic interaction with a weak H-bond (w-w pairing) explains the inhibitory action of antagonists on factor Xa activity ³⁷⁰. The compound with a quaternary ammonium cation [–N(Me)₃⁺] is ~1100-fold more active than the ammonium ion (–NH₃⁺) compound, indicating that hydrophobic pairings with weak (w-w) H-bonds are more favorable than the s-w interactions. Although the –NH₃⁺ may not interact with the aromatic rings, logic analysis indicates that –NH₃⁺-π interactions are less favorable than –N(Me)₃⁺-π interactions. Also, van der Waals interactions are insufficient to explain the large potency differences. Thus, mixed s-w H-bond pairings can decrease protein-ligand binding affinity even when H-bonds are strong, offering new mechanistic insight into why some strong H-bonds do not enhance ligand binding affinity; generalizations that H-bonds contribute minimally to binding affinity may therefore be inaccurate. Notably, s-w H-bond pairings with a ligand atom can significantly reduce its binding affinity as is demonstrated by the relatively weak binding affinity of the scytalone dehydratase inhibitor 62 (Figure 16) ³⁶⁰, with its apolar hydrogen atom (H_L = 0 < H_W) interacting with two strong H-bond donors from Tyr30 and Tyr50 (H_P = 21.6 > H_W). The s-w pairing interaction decreases binding affinity by ~90-fold.

7. Prediction of Hydrogen Bonds Formation

The prediction of hydrogen-bond (H-bond) acceptor ability is crucial in drug design. This important property is quantified by Graton et al. ³⁷¹ in a large pK_BHX database of consistently measured values. Their aim was to expand the chemical diversity of the studied H-bond acceptors and to increase the range of H-bond strength considered. Two quantum chemical descriptors were contrasted, called ΔE(H) (the change in the energy of a topological hydrogen atom upon complexation) and V_min (the local minimum in the electrostatic potential on the H-bond accepting site).

Hydrogen bonding ^{128, 205, 345, 372, 373, 374} is a key element of molecular recognition ^{368, 375} and is implicated in diverse physicochemical phenomena such as crystal packing ^{376, 377, 378}, DNA base-pairing ^{379, 380}, protein folding ^{1, 381, 382, 383, 384}, and enzyme specificity ³⁸⁵. The cohesiveness of liquid water that drives hydrophobic association ^{386, 387} in aqueous media is a consequence of cooperative hydrogen bonds between water molecules. In the pharmaceutical context, hydrogen bonding influences solubility of drugs, both in water and in lipids, and the affinity with which they associate with their targets. Partition coefficients ^{388, 389}, especially when measured for the 1-octanol/water system, are widely used in drug discovery, and differences between values measured for a compound in different systems (e.g., cyclohexane/water) reflect the hydrogen-bonding characteristics of the compound ^{241, 358, 375, 390, 391, 392, 393}. Hydrogen bond (HB) strength can be quantified by the stability of the 1:1 complex between an HB donor and an HB acceptor in a nonpolar solvent such as tetrachloromethane or 1,1,1-trichloroethane ^{25, 76, 394, 395}. The 1:1 complex stability measured using a reference HB donor such as 4-fluorophenol ^{76, 394, 395} or 4-nitrophenol ²⁵ can be used to define scales for HB basicity, and a reference HB acceptor such as 1-methylpyrrolidone can be used to quantify HB acidity in an analogous manner ²⁵. Measured HB basicity (and acidity) can be used to guide molecular design, although there are limitations to how the measurements can be used. First, measurements made for compounds which have nonequivalent HB acceptors cannot, in general, be used to characterize the different HB acceptor sites for a molecular structure, although in some cases ³⁹⁶ the individual isomeric 1:1 complexes can be observed. Second, in polar solvents such as water or in ligand-protein complexes, formation of multiple interactions can perturb HB basicity and acidity of individual acceptors and donors. This is a particular issue for carbonyl oxygen atoms, which are typically associated with two HB acceptor sites, and hydroxyl groups for which an HB donor and acceptor are in close proximity ^{358, 397}. Outside narrowly defined structural series, pKa is not an effective predictor of HB acidity or basicity ^{25, 395}. Modeling HB Basicity with V_min. Many medicinal chemists will be unfamiliar with the MEP minima used to model HB basicity so a good way to start this section is to note that these can be considered to be roughly equivalent to “lone pairs”. The Vmin value at the MEP minimum can be thought of as the “strength”, “intensity”, or “availability” of the “lone pair”.

V_min is a more direct (and less arbitrary) measure than an atomic charge of what an HB donor would “see”. MEP minima typically lie within the van der Waals molecular surface, and computational chemists should be aware that these minima cannot, in general, be reproduced using atom-centered multipoles. When only a single MEP minimum is present in the molecular structure, HB basicity can be modeled ³⁹⁸ according to a linear model with intercept A and slope B:

(24)

If two HB acceptor sites are present in the molecular structure, the observed association constant K is simply the sum of the association constants (K₁ and K₂) for the individual acceptors ^{395, 398}:

(25)

There are two simplifying assumptions that can sometimes be made in order to fit measured pK_BHX when two HB acceptor sites are present in a molecular structure. First, when the HB acceptors are equivalent (K₁ = K₂), the measured pKB_HX can be corrected statistically and used for fitting as if only a single HB acceptor is present. Second, when one of the HB acceptors is likely to be much stronger than the other (e.g., oxazole), we can neglect the contribution of the weaker acceptor to HB basicity when fitting pK_BHX to V_min. In this scenario, the relevant V_min values are useful for testing the validity of the assumptions.

Oxygen, sulfur, and fluorine atoms are typically associated with two or three MEP minima with comparable V_min values that are not in general equivalent, and it is not, in general, valid to statistically correct pK_BHX in these situations. One solution to the problem of fitting measured pK_BHX to comparable, nonequivalent V_min values is to use nonlinear regression as exemplified by eq 26 in which V₁ and V₂ are two nonequivalent V_min values:

(26)

Although eqs. 25 and 26 can be generalized to any number of HB acceptor sites, the pK_BHX value measured for a compound with multiple, nonequivalent HB acceptors has no thermodynamic significance ³⁹⁵. Furthermore, the sum of association constants for all isomeric 1:1 complexes is not directly relevant to the solvated state in which all HB acceptors and donors of the solute can interact simultaneously with solvent molecules.

HB basicity measurements for compounds with nonequivalent acceptors are only useful if the association constants for individual acceptors can be quantified ³⁹⁶, and this is not generally the case for compounds of interest to medicinal chemists. Predictive models usually represent the only practical means of quantifying HB basicity of individual acceptors in molecular structures of interest to medicinal chemists, although measurements made for structural prototypes are absolutely essential for developing the predictive models.

Shalaeva et al. study ²⁴¹ demonstrates that ΔlogP_oct-tol (difference between logP_octanol and logP_toluene) describes compounds propensity to form intramolecular hydrogen bonds (IMHB) and may be considered a privileged molecular descriptor for use in drug discovery and for prediction of IMHB in drug candidates. They also identified experimental protocols for acquiring reliable ΔlogP_oct-tol values on a set of compounds representing IMHB motifs most prevalent in Medicinal Chemistry, mainly molecules capable of forming 6-, 7-member IMHB rings. And furthermore, computational ΔlogP_oct-tol values obtained with COSMO-RS software provided a good estimate of experimental results and can be used prospectively to assess IMHB.

The incorporation of an intramolecular hydrogen bond (IMHB) into a molecule is gaining a great deal of interest in drug design as indicated by the number of papers recently published in key Medicinal Chemistry journals. ^{209, 217, 247, 399, 400} The presence of IMHB has been shown to significantly alter molecular properties due to formation of various conformers that in turn influence solubility, permeability, PK/ PD processes, and protein binding affinity ^{207, 210, 401, 402}.

The IMHB as described by Desiraju ¹²⁴ is an attractive interaction in which an electropositive hydrogen atom intercedes between two electronegative fragments of the same molecule and holds them together. A hydrogen bond is strong enough to restrict rotation of fragments by forming most commonly ^{210, 217, 401, 402} membered rings. Importantly, IMHBs are weak enough to allow these fragments to come apart and lose their orientational specificity in high dielectric media such as water. The chameleon like nature of an IMHB becomes apparent when one realizes that in water an IMHB is unlikely to form and the polar groups may serve to increase solubility by readily forming intermolecular hydrogen bonds with water. Alternatively, molecules that can participate in IMHB shed water more readily when entering a low dielectric environment like a hydrophobic phospholipid bilayer. In this circumstance IMHB results in lipophilic, less polar molecular conformations which are expected to have higher passive membrane permeability ⁴⁰³. In other words, a decrease in polarity is sometimes achieved through the formation of IMHBs, where the hydrogen bond donor (HBD) and acceptor (HBA) atoms are effectively shielded from water, thereby reducing the energetic penalty of desolvation required in moving from an aqueous environment through a phospholipid bilayer ⁴⁰⁰.

The consequences of IMHBs to medicinal chemists are significant but often under-recognized and seldom predicted. For instance, lipophilicity may be underestimated when determined by calculated logP (clogP) in molecules with IMHBs, while hydrogen bond donor and acceptor counts are overestimated. Additionally, clogP, as well as hydrogen bond donor and acceptor counts, are part of the ubiquitous RO5 parameters ¹⁹³, used to predict drug like properties and permeability. When IMHB are present these RO5 counts can be effectively stretched, broadening drug like property space allowing more diverse drug design ^{207, 404}. Likewise preferred property space for Central Nervous System (CNS) drugs may be extended when IMHB are present, as hydrogen bond donor count and clogP are both parameters in the CNS Multi-Parameter Optimization (CNS MPO) score ⁴⁰⁵. In support of this notion it was also found that ΔlogP_oct-alk correlates with brain penetration and oral absorption ^{406, 407}.

Unfortunately, one cannot simply examine a given 2D structure and immediately delineate the presence of one or more IMHBs and determine their strength because the thermodynamic equilibrium of closed versus open conformations depend on a number of complex factors (e.g. geometry, type of solvent and others) ⁴⁰⁸ acting simultaneously. The most common tools used to investigate IMHBs are spectroscopy (NMR, infrared and Raman, microwave), diffraction (X-ray and neutron diffraction), calorimetry and theoretical methods ⁴⁰⁹. However, many of these techniques are not high throughput and data produced often require detailed interpretation by experts. These issues lead us to look for additional methods.

LogP is one of the most widely used parameters in drug design and it has been considered for evaluation of IMHB ^{391, 402, 410}. It has been demonstrated using solvatochromic equations that the difference between logP values obtained in different biphasic systems (ΔlogP), for example octanol/water and alkane/water (ΔlogP_oct-alk = logP_oct – logP_alk), is informative of IMHB when the solvents are very different from each other ³⁵⁶. The idea that ΔlogP_oct-alk is informative of IMHB and the reports that ΔlogP_oct-alk correlates with brain penetration and oral absorption ^{406, 407} lead to the exploration of ΔlogP_oct-tol (ΔlogP_oct-tol = logP_oct – logP_tol ^{358, 411, 412, 413, 414}. The investigation of IMHB by ΔlogP was proposed some time ago ⁴¹⁵, however this approach was not widely implemented mainly because the practical tools, both experimental and theoretical, to obtain logP_alk data for large series of compounds were limited.

The main goal of this study is to demonstrate that ΔlogP_oct-tol (ΔlogP_oct-tol = logP_oct – logP_tol ) distinguishes compounds with high propensity to form IMHB and to develop a protocol for its implementation in active Medicinal Chemistry projects where series of similar compounds are often available for relative comparisons.

According to the results reported by Shaaleva et al. ²⁴¹, ΔlogP_oct-tol should be included in Medicinal Chemistry design for optimization of physical chemical properties, as a privileged molecular descriptor for delineating the propensity of compounds to form IMHB.

No general guidelines are reported in the literature on the interpretation of ΔlogP in relation to the presence and the strength of IMHBs. Therefore the analysis of ΔlogP data was aimed at obtaining an IMHB interpretation method.

As discussed earlier, it is assumed that toluene, similar to apolar solvents, promotes folded conformations and formation of IMHB when possible, whereas the reverse is true for molecules in water and, to a lesser degree, in octanol. Therefore, the difference between logP_oct and logP_tol (i.e. ΔlogP) should reflect the propensity of a compound to form IMHB. The ΔlogP value by itself does not indicate the formation of IMHB.

8. The Importance in Drug Design of Similar to HB Bonds

As previously commented, non-covalent interactions, especially hydrogen bonds (HBs), are known to be responsible for the conformation and 3D structure of biomolecules like proteins and DNA. Among the so-called non-conventional hydrogen bonds ⁴¹⁶, there is the R-C-H---X group (X=O, Cl and N). These are weak hydrogen bonds ^{17, 373}, since the energy of these secondary bonds is smaller due to the lower acidity of the involved hydrogen atoms. These non-covalent interactions as halogen bonds ^{2, 3, 4}, pnicogen bonds ^{5, 6} and tetrel bonds ^{7, 8} can contribute as well to the stability of certain molecular conformations. The study of intramolecular interactions is very important in the design of pharmaceutical drugs, particularly in the context of conformationally flexible molecules. Conformation-controlling intramolecular interactions in drug molecules have a direct influence on the binding modes of the drugs with the respective targets ^{53, 54, 55}. In particular, intramolecular peri-interactions have been widely studied in the literature in naphthalene and other related systems ^{56, 57, 58, 59, 60}.

Chalcogen bonds ^{417, 418, 419, 420, 422} (YB) are one of the less studied non-covalent interactions. However, a few articles have been devoted in the literature to the study of intramolecular YBs. The nature of the YB has been rationalized by Politzer et al. based on the σ–hole concept ^{423, 424, 425}. The term σ –hole refers to the electron-deficient outer lobe of a p orbital involved in forming a covalent bond. Therefore, in YBs, as well as in other non-covalent interactions, the importance of the electrostatic interaction term is uppermost ^{426, 437}. Several works have been carried out regarding the tunability of the mentioned σ –holes in different types of interactions, including halogen, chalcogen, and pnicogen ^{438, 439, 440, 441}. In particular, the effect on the σ –hole upon substitution on aromatic rings has been previously studied ^{438, 442, 443, 444}.

The relationships between the interaction strengths and various geometric and electronic structures, as well as topological properties, have been established ⁴⁴⁵, with a particular focus on the uniformity of hydrogen bonds and halogen bonds types of interaction. The dependence of the BSSE-corrected interaction energy (ΔE_cor) on the interatomic distance (rA_/N) appeared to be nonlinear for both halogen-bonded and hydrogen-bonded systems; the relationship between ΔE_cor and the difference between r_A/N and the sum of the van der Waals radii (Δr_A/N) can be fitted to a combined quadratic regression equation. Furthermore, it has been demonstrated that the linear correlations between ΔE_cor and r_b(BCP) (the electron density at bond critical points in the A/N bond) and its Laplacian V²r_b(BCP) can be used to provide a combined description of hydrogen bonds and halogen bonds ⁴⁴⁵.

Strong halogen-bonded complexes were found to exhibit different linear correlations from weak halogen-bonded and hydrogen-bonded systems ⁴⁴⁵. Nevertheless, for the latter two types of system, a uniform regression equation can be constructed. These relationships not only improve our understanding of the nature of halogen bonding but also provide a feasible approach for predicting or determining the relative strengths of hydrogen bonds and halogen bonds, in particular when both types of non-covalent interaction coexist and compete with each other.

A theoretical study of the peri interactions, both intramolecular hydrogen (HB) and chalcogen bonds (YB), in 1-hydroxy-8YH-naphthalene, 1,4-dihydroxy-5,8-di-YH-naphthalene, and 1,5-dihydroxy-4,8-di-YH-naphthalene, with Y = O, S, and Se was recently carried out ⁴⁴⁶, stating that the systems with a OH:Y hydrogen bond are the most stable ones followed by those with a chalcogen O:Y interaction, those with a YH:O hydrogen bond (Y = S and Se) being the least stable ones. The electron density values at the hydrogen bond critical points indicate that they have partial covalent character. Natural Bond Orbital (NBO) analysis shows stabilization due to the charge transfer between lone pair orbitals towards empty Y–H that correlate with the interatomic distances. The electron density shift maps and non-covalent indexes in the different systems are consistent with the relative strength of the interactions. The structures found on the CSD were used to compare the experimental and calculated results.

In a study calculating the equilibrium conformer compositions for 2X-ethanol and 2X-phenol (X = F, Cl) in solution ¹⁰⁴, the predicted OCCF gauche/trans ratio for 2F-ethanol was well reproduced in comparison with available experimental compositions. The predominant gauche structure maintains an intramolecular hydrogen bond in carbon tetrachloride (HB structure), whereas HB and NoHB gauche conformers appear in nearly the same fraction in aqueous solution. The internally hydrogen-bonded conformer is predominant also for 2X-phenol species. Solution structure modeling predicts weak hydrogen-bond formation capacity for both the covalently bound F and Cl atoms, even in conformations where they are fully exposed to hydration.

9. Conclusions

In this review, the effect of weak intermolecular interactions on the binding affinity between ligand-protein complexes in order to improve drug efficacy has been explored with special emphasis in the important features of the role of hydrogen bonding in the drugs activity. Studies of hydrogen bonds were presented and discussed as they raise more questions than have been previously answered. For example, if energetically stabilized drug-like compounds are trapped at the bioactive core of the target site, holding all the biochemical and conformational features, then how is the side effect manifested? How is the binding affinity between drug-ligand complexes associated with drug efficacy? If a suitable environment is provided, are hydrophobic interactions and hydrogen bonds interchangeable? Furthermore, computational biologists (in-silico) are challenged to find supporting factors that bring long-range associated ligand-target complex molecules into small regions where biological activity can be altered. To resolve all these issues, a multi-model approach is needed that explores the dynamic nature of weak intermolecular interactions at the target-drug interface.

Almost all modern computational chemistry techniques are based on the analysis of molecular properties, but very few address the interactions between different parts of the molecule and the interactions of molecules with one another. In this review, we have demonstrated that rules characterizing H-bonding in drug molecules and its targets can be induced from properly preprocessed 2D structural data, because it has been clearly demonstrated in this study that prediction of hydrogen bond energies is possible.

The fundamental function of the intramolecular hydrogen bonding in drugs has been highlighted in this context as well as the differences arising in the interaction ligand-protein when tautomerism exists and hydrogen bonding patters may be altered, as well as the binding affinity between ligand and protein.

Inspection of the outliers in this study proved very stimulating. For instance, it provides pattern recognition opportunities and molecular design inspiration to medicinal chemists and machine learning methods. Perhaps more importantly, it is a sober reminder of the impact that minimal chemical modification can have on biological activity, even when taken outside of the IMHB context discussed here.

Competing Interests

The author declare no competing financial interest.

Acknowledgements

The autor thanks the MCINN Project Consolider Ingenio CSD2010-00065 for financial support.

References