We first evaluated by HINT the free energy of binding for 53 protein–ligand complexes formed by 17 proteins of known three-dimensional structure and characterized by different active site polarities (Marabotti et al., 2000; Cozzini et al., 2002). This analysis was carried out without the contribution of bound water molecules to the binding energy. A successive analysis considered such contributions (see Section 4.1) (Cozzini et al., 2004; Fornabaio et al., 2004; Kellogg et al., 2004). To demonstrate that HINT was able to correctly predict the interaction energy between protein and ligands, independently of ligand and protein active site polarity, the selected protein’s active sites vary from very apolar, such as retinol-binding proteins, where the hydrophobic contribution and, consequently, the entropic contribution to the binding energy is predominant, to polar, such as penicillopepsin, where the free energy of interaction is mainly associated with contributions from Coulombic interaction between polar or ionizable residues. Another key feature of the analyzed set of protein–ligand complexes is represented by the large range of binding strength that varies over nine orders of magnitude, with ΔG ranging from −2 to −15 kcal/mol, as calculated by either inhibition constants or dissociation constants (Cozzini et al., 2002).

HINT scores for the protein–ligand complexes were plotted against the experimental binding affinity, obtaining a remarkably good correlation with a standard error of 2.6 kcal/mol, which translates into a prediction of affinity within about two orders of magnitude (Figure 2). Even better standard errors (1 kcal/mol) were obtained within single sets of specific ligand–protein complexes, such as trypsin, thrombin, and tryptophan synthase. Given the polarity heterogeneity of the 53 protein–ligand complexes, the prediction was very good.

A key requirement for obtaining good correlation between computational and experimental data, thus predicting binding affinity for ligands with known three-dimensional complex structures, is the homogeneity between the pH at which binding affinities are measured in solution and the pH at which structures have been determined. In fact, frequently, solution data are determined at pH values quite different from crystal structure medium. As a result, the protonation state of groups might be incorrectly attributed (or modeled), thus affecting the scores. A telling story is represented by an analysis of the binding affinity of penicillopepsin–ligand complexes measured at pHs 3.5, 4.5, and 5.5 and by three-dimensional structures determined at the same pH values. HINT score prediction using these data produced a linear correlation with r ² equal to 0.99 (Cozzini et al., 2002), indicating that when pH values of solution and structure data match each other, the HINT prediction is remarkably good.

To address, in a more general way, the dependence of the HINT score on the protonation state of interacting groups of ligand and protein, a protocol called “computational titration” was designed (Fornabaio et al., 2003; Kellogg et al., 2004; 2006; Spyrakis et al., 2004). The steps required by this method (Figure 3) are as follows: i) the identification by a careful inspection of the three-dimensional structures of the ionizable groups of ligand and protein that can contribute to the free energy of binding; ii) the generation of all possible ionization state models for the interacting groups; iii) the energy minimization of the generated models with no movement of the heavy atoms (isocrystallographic models); and iv) the determination of the score for each model with HINT. This is a rather simple problem when only a few ionizable residues or functional groups are involved but increases rapidly in complexity. Each acid has three possible states (unprotonated and protonation at the two oxygens), while each amine or thiol has two, and—if considered—guanidyl could have five. The result is the dependence of the interaction free energy on group ionization states and the identification of the pattern of group ionization states, which leads to the optimal interaction strength. This, in turn, can be correlated with a well-defined pH value (Spyrakis et al., 2004).

We applied the computational titration algorithm to the analysis of the interaction between neuraminidase and nine inhibitors for which three-dimensional structures and inhibition constants were known from the literature (Fornabaio et al., 2003), between HIV protease and its substrate peptide (Spyrakis et al., 2004) and between dihydrofolate reductase and three ligands (Kellogg et al., 2006). In the case of the computational titration procedure for neuraminidase, one proton at a time was introduced into the molecular model. Here, the protein possesses three ionizable amino acid sidechains, namely, Glu 119, Asp 151, and Glu 276, in the active site and one carboxylic moiety on one of the ligands. Thus, a total of four protons were added to the models, plus the starting 0-level model in which all ionizable groups are unprotonated. For each added proton, all possible models were generated, some of them being more chemically probable than others. As shown in Fornabaio et al. (2003), the dependence of HINT score vs. added protons exhibits, not surprisingly, a bell shape, suggesting the optimal protonation state of the ionizable residues, and thus the pH for optimal interaction strength. A further step in the optimization of the computational titration, which requires building and energy minimizing sometimes thousands of models with specific protonation for each protein–ligand complex, was to exclude high-energy models (chemically unplausible) and include a statistical mechanics evaluation of all models. This refined approach was applied to the computational titration of HIV protease and its peptidic substrate (Spyrakis et al., 2004). There are four ionizable residues in the active site, i.e., Asp 25α, 29α, 30α, and 25β, and four potentially interacting ionizable groups on the substrate, i.e., three carboxylates and one amine. In this case, 4,374 unique protonation models were generated, ranging from the most basic with all sites deprotonated (overall charge of −7) to the most acidic with all sites protonated (overall charge of +1). Titration analysis showed: i) a quite large range of HINT scores for models containing the same number of protons indicating favorable and unfavorable proton distributions (Spyrakis et al., 2004), ii) a sharp increase in the HINT score as protons were added that levels off upon five added protons, and iii) a good correlation between experimental data, collected between pH 5 and 3, and the Boltzmann-averaged HINT scores.

Another factor that shows a profound impact on the predictions based on HINT, but also to all other codes, is the quality of the three-dimensional structures. The higher the resolution, the lower the standard error in the HINT prediction. We have observed that by selecting within the 53 complexes only structures determined with a resolution better than 2.5 Å, the standard error of the prediction improved to 1.8 kcal/mol compared to 2.6 kcal/mol when the included complexes possess a resolution within 3.2 Å (unpublished data). A better crystallographic resolution leads to a more precise geometry and well-defined distances of interacting groups that impact HINT scores.

Overall, the HINT score for a protein–ligand complex depends, in addition to other contributions, on hydrogen-bonding contribution. A particular property of HINT is that it is very sensitive to the positioning and orientation of hydrogen-bonding protons. It is well known that hydrogen atoms can rarely be detected crystallographically, except for structures with resolutions higher than ∼1 Å. Thus, automated procedures were put in place to insert hydrogen atoms bound to heavy atoms and to correctly orient them to optimize the geometry for hydrogen bonds and, therefore, their strength and scores, without altering the position of the crystallographically determined heavy atoms. The energetic contribution of hydrogen bonds also plays a significant role when water molecules bound to protein active sites and/or ligands are considered (see the following section).

The literature provides numerous cases in which water plays crucial roles in mediating protein–ligand interactions. Examples are represented by bosutinib binding to Src kinase (Levinson and Boxer, 2014), thermolysin and its water network stabilized by carboxybenzyl-Gly-(PO₂)-L-Leu-NH₂−based ligands (Krimmer et al., 2014), the water networks of the adenosine A2A receptor and their perturbations resulting from ligand binding (Congreve et al., 2012), or by the well-known case of HIV-1 protease (Lam et al., 1994). Trying to rationalize the associated water roles, we used HINT to estimate the energetics of binding site waters in 23 HIV-1 protease/inhibitor complexes that belong to the different classes of hydroxyethylenic ligands, peptidomimetic diol derivatives, cyclic ureidic derivatives, and cyclic sulfamide derivatives. The latter pair was specifically designed to displace a conserved water molecule (water 301) in the binding site. A careful analysis of the X-ray structures of the complexes showed the presence of three different structural water categories, that is, water 301, symmetrically located in the binding site with respect to the catalytic one and H-bonding Ile 50 and Ile 150 backbones, waters 313 and 313′, located in a more peripheral region, and waters 313bis and 313bis’, inserted in small clefts in the binding site. The protein–ligand interaction calculated by HINT, with the inclusion of the water contribution, showed that bridging water molecules, such as 301, significantly contributes to the overall energy (about 4–6 kcal/mol) and leads to a better correlation between experimental and predicted free energy. In contrast, the contribution of more peripheral waters is variable or not significant and strongly dependent on the chemical nature and size of the ligand (Fornabaio et al., 2004). However, the binding energy gained by the new ligand-to-protein (Ile 50, Ile 150) interactions in the water 301-displacing cyclic derivatives was very similar to the energy previously associated with the water-to-ligand and water-to-protein interactions (Figure 4). As a conclusion, we recommended including the contribution of water molecules when predicting the free energy of binding but paying specific attention to ensure that this adds information and not just noise.

We implemented the previous work by evaluating not only the energy associated with protein–water interaction but also the geometry quality of the H-bonds formed by each water molecule to obtain a more reliable estimation and prediction of water’s role in mediating protein–ligand association (Amadasi et al., 2006). First, 12 solvated proteins, in both the apo and holo form complexed with inhibitors, were considered for a total of 2,186 co-crystallized water molecules. Water molecules were classified and separated into the following categories: 1) water molecules in active sites; 2) water molecules deeply inserted into cavities; 3) buried water molecules; 4) first shell external water molecules (within 4 Å of the nearest protein heavy atom); and 5) second shell external water molecules (more than 4 Å from the protein; see the study of Amadasi et al. (2006). Each water molecule was evaluated and classified according to the interactions made with the protein, by means of the HINT score, and to the number of potential H-bonds, by means of the rank algorithm (Kellogg and Chen, 2004).

The latter can evaluate the count, strength, and geometry of potential hydrogen bonds for each water molecule in the cognate structure and is also used to optimize water position and hydrogen placement. Rank can vary from 0, for waters not involved in any hydrogen bonds, to around 6, for waters forming four high-quality hydrogen bonds in terms of bond length and angle geometry, and it is calculated by the following equation: R a n k = ∑ n 2.80   Å / r n + ∑ m cos Θ T d − Θ n m / 6 , (3) where r _n is the distance between the water oxygen and the target heavy atom; Θ _Td is the ideal angle of 109.5°, and Θ _nm is the angle between targets. Also, any angle less than 60° is rejected, as well as the corresponding bond.

The analysis showed that both HINT and Rank scores increase when evaluating water molecules in the second and first hydration layers, to waters in active sites, up to waters in protein cavities, or completely buried in the protein matrix.

Then, 15 protein–ligand complex binding sites, in which the presence of at least one bridging water molecule was reported in the literature, were analyzed to provide reference HINT scores and Rank values for bridging waters. We observed that the average Rank for protein–water and ligand–water was estimated to be 3.0 and 1.5, respectively, suggesting that proteins are better able to embed bridging waters than ligands, which is quite reasonable considering the residue sidechain flexibility and the presence of clefts and dips in active sites. We have to consider that a Rank equal to 3 might correspond to three H-bonds but also to two very good bonds in terms of distance and geometry. However, no significant difference was provided in terms of the HINT score for protein–water and ligand–water interactions.

Finally, we analyzed a set of nine proteins in both native and complexed state and classified water molecules in active sites in the following categories according to their relevance: 1) conserved water molecules in binding sites bridging protein–ligand interaction; 2) conserved water molecules in binding sites not relevant for protein–ligand interaction; 3) conserved water molecules in binding site cavities; 4) conserved water molecules in peripheral binding site regions; 5) water molecules displaced by ligands replacing their function, i.e., functionally displaced; 6) water molecules displaced by ligands only occupying their room, i.e., sterically displaced; and 7) missing waters. Water molecules with Rank >1.5 and HINT score <150 can be considered as sterically displaceable and, thus, potentially removable by ligands at moderate cost. Water molecules having, instead, Rank values in the 1.5–4.0 range and HINT scores >150 can be more relevant and quite likely retained with respect to drug design, while water molecules with Rank >4 are too buried in the binding site and of less interest (Amadasi et al., 2006).

The diagnostic potential of the HINT score and Rank in predicting displaceable, bridging, or buried water molecules was further implemented in the more quantitative Relevance metric (Amadasi et al., 2008). Although water molecules are usually excluded from consideration in structure-based drug design, the rapid identification of Relevant water molecules that might be retained or displaced by ad hoc designed ligands (any water can be, in principle, displaced) can significantly affect the energetics of protein–ligand interaction. By using a training set of thirteen proteins, crystallized in the native and complexed state, with a total number of 125 water molecules, we combined information coming from the HINT score and Rank in the following equation, describing the overall probability (P _A ) of specific water to be Relevant: P A = P R W R + 1 2 + P H W H + 1 2 W R + 1 2 + W H + 1 2 , (4) where P _R and P _H are calculated by nonlinear polynomial regressions and correspond to the percent probability for conservation based on the Rank and HINT score, respectively, while W _R and W _H are the weights of Rank and HINT score probabilities, respectively (Amadasi et al., 2008). This analysis showed Rank and HINT score values associated with conserved waters: P _R > 60%, corresponding to Rank >2.3, and P _H > 80%, corresponding to HINT score >400. Otherwise, P _H > 40%, although corresponding to HINT score >100, was considered indicative of non-conservation. When applied on a test set of nine native and complexed proteins, for which specific bridging water molecules were previously identified, 59 of 68 water molecules were correctly predicted, corresponding to a success rate of 87%, which increased to 92% when only water molecules from X-ray structures having a resolution <2.0 Å were considered. Importantly, the crystallographic B factors, often used in such analyses did not improve this model.

In a later contribution, we designated water molecules in protein matrices as “cold” and “hot” according to their internal energy and possible displacement (Spyrakis et al., 2017). Hot water can also be considered “unhappy” water, that is, not stable in binding sites or at protein surfaces, because of the presence of extensive hydrophobic regions. Being unhappy, they would benefit from being displaced toward a more polar environment as the bulk, where they could maximize the number of H-bonds and provide a favorable variation of the binding free energy. Cold water, being stably bound in polar environments can be, instead, considered as “happy” molecules able to extend protein/ligand functions and to participate in the protein structure, dynamics, and function. Their displacement might not be at a trivial expense and may be associated with meaningless or even unfavorable variation of the binding free energy.

Behind the hot/cold classification is, in fact, their enthalpic and entropic contributions. Several considerations can be drawn that might be useful with respect to a drug design perspective (Ahmed et al., 2013; Spyrakis et al., 2017): i) hot water molecules in binding sites can be easily displaced with a manifest gain in entropy and a negligible loss in enthalpy; ii) attention must be paid when displacing cold water molecules since the entropic gain might not balance the enthalpic loss; iii) the enthalpy/entropy compensation associated with cold water displacement can lead to quite small changes in the binding free energy, thus often making drug design efforts ineffectual; iv) hot water networks around protein–ligand complexes hide hydrophobic moieties stabilizing the association; and v) both hot and cold water molecules at the protein–protein interface can provide precious indications when designing protein–protein inhibitors.

Clearly, the emphasis on modern drug discovery has shifted over the past decade or so toward higher-throughput modeling approaches, such as virtual screening. Before using HINT for extensive docking experiments, we were curious whether docking scores from various scoring functions correlated better with RMSD (root-mean-squared distance) or free energy of binding. In other words, does reproducing a crystal structure by docking, which is how most scoring functions are optimized, necessarily produce accurate predictions of binding energy? We examined 19 protein–ligand complexes for which X-ray crystallographic structures and binding energy data were available, and we calculated experimental vs. computationally-derived free energy correlation by means of the HINT free energy scoring tool and other scoring functions (Spyrakis et al., 2007a). Correlations drawn after re-docking the cognate ligands in their corresponding protein structures were generally better with the HINT scoring function (Spyrakis et al., 2007a). Also, it was obvious from this study that scoring functions uniquely calibrated for the dataset or sets under study should nearly always be preferable to universal scoring functions.

The HINT scoring function described earlier has indeed been a powerful tool for virtual screening (Salsi et al., 2010; Farzan et al., 2011; Spyrakis et al., 2014; Obaidullah et al., 2018; Kayastha et al., 2022) with a few caveats, including: 1) the scoring function is extremely sensitive to structure, especially properly optimized hydrogen-bonding interactions that include the donor hydrogen’s position, and 2) while it is more than fast enough for scoring docking models, it is not, at present, in the “giga” docking timescale range. Thus, HINT, in our hands, has been used as a second scoring function after docking with programs like GOLD or AutoDock and their native scoring functions have been applied as primary filters. In this use case, our in-house studies generally have been relatively small-scale and speculative, i.e., purchasing 10–100 compounds for actual assay with a variety of low-to-medium-throughput methods. Nevertheless, about one-in-five of the compounds identified by virtual screening with HINT scoring were active.

Three examples of this research are as follows: 1) 14 pentapeptides, binding to Haemophilus influenzae O-acetylserine sulfhydrylase (OASS) with affinities ranging from μM to mM, were identified with good correlation between computational and experimental data, excluding peptides bearing a positive charge, which are likely overestimated by HINT. These results, combined with X-ray structures of the three best complexes, defined a pharmacophoric scaffold for the design of peptidomimetic inhibitors for this enzyme (Salsi et al., 2010). 2) To identify new antimicrobials toward Treponema denticola cystalysin, we performed virtual screening on 9,357 compounds with the FLAPsite algorithm (Baroni et al., 2007), docked those most promising, and rescored with HINT in a consensus approach. Among the 17 compounds selected for testing, two showed IC₅₀s in the low μM range, identifying interesting hits for the development of novel antimicrobials (Spyrakis et al., 2014). 3) The most recent publication using HINT in this way (Kayastha et al., 2022) details the discovery of three new potent eIF4A inhibitors that diminished the viability of diffuse large B-cell lymphoma (DLBCL) cells. The actives were discovered via target-based virtual screening in the MolPort database, followed by pharmacophore-based and analog screening. Modeling suggests that these compounds clamp eIF4A into mRNA in an ATP-independent manner and depress eIF4A-dependent oncogene expression. Eukaryotic translation initiation factor 4A is essential in translation initiation as it unwinds the secondary structure of messenger RNA upstream of the start codon and enables active downstream ribosomal recruitment.

While protein–ligand interactions are complex on their own, the addition of more degrees of freedom when two proteins interact is another level of complexity. For example, in a typical protein–ligand interaction, there will likely be only a few relevant water molecules, but at a protein–protein interaction surface, there can be one or two dozen. Using the aforementioned computational titration algorithm (Fornabaio et al., 2003; Spyrakis et al., 2004; Kellogg et al., 2006) would typically involve 2–3 ionizable residues or functional groups in a protein–ligand complex, but again, at a protein–protein surface, this could be a much larger number.

To explore the water issue, we adapted the rank and relevance algorithms to protein–protein interactions in order to probe the roles of water molecules found specifically at protein–protein interfaces (Ahmed et al., 2011; 2013). The result was the surprising conclusion that less than one-quarter of waters found at such interfaces were engaged in truly structure-supporting bridging roles, and just over one-quarter of them had unfavorable interactions with respect to both proteins. The remainder, over one-half, had favorable interactions with one protein and unfavorable interactions with the other (Ahmed et al., 2011). We proposed a new structural motif—the hydrophobic bubble—to describe the phenomenon of one or a collection of water molecules “trapped” in a hydrophobic cavity formed by the protein–protein association (Figure 5). Moreover, we hypothesized that these bubbles are critical to protein function as points/regions of instability to assist in protein–protein dissociation when needed. Two implications of these observations are worth mentioning: 1) it will be very difficult to predict the presence of water in such loci based only on energetics as they are by definition unstable, although water clusters do gain enthalpic advantage through their interactions with each other, and 2) these phenomena may present a unique opportunity for drug discovery at protein–protein interfaces as they are the true “hotspots” (Spyrakis et al., 2017).

To further examine our insistence on the importance of water at protein–protein interfaces, we reported a docking study on a small number of protein–protein complexes both “dry”, as allowed by the native ZDOCK (Pierce and Weng, 2008; Pierce et al., 2011) program, and “wet”, where we tricked ZDOCK into recognizing interfacial water molecules (Parikh and Kellogg, 2014) (Figure 6). The latter models were far superior by HINT score and numerous metrics developed by the Critical Assessment of PRedicted Interactions (CAPRI) communitywide experiment (Janin et al., 2003; Mendez et al., 2005), but creating such models was very tedious at that time due to the limitations inherent in the ZDOCK code, e.g., its automatic deletion of anything that was not in its library of 20 amino acids. This, coupled with the inherent problems in predicting the presence of water molecules that are “not” energetically viable, although obviously present, set us in a different direction with respect to HINT-based structure prediction (section 5.3), including protein–protein docking.

As we were developing HINT, implementing an intramolecular score was a very simple addition. However, it was not immediately obvious how to apply it to interesting and important problems until we started thinking about the difficulties of dealing with low-resolution crystal structures and how their reflection data are collected and refined. When the resolution is high, the reflection data are more than adequate to model the structure within the electron density envelopes and features, such as hydrophobic interactions, are seen. However, at low resolution, somewhat crude molecular mechanics forcefield algorithms are applied, and hydrophobic interactions—that are not explicit in these force fields—are often lost. We used the intramolecular HINT score as an adjuvant to contemporary structure refinement protocols (Koparde et al., 2011). We showed that low-resolution structures refined by including our HINT intramolecular score-based protocol were significantly more native-like based on structure quality metrics.

More recently, Agosta et al. (2022) used the intramolecular HINT score to evaluate the stability of proteins in response to single (and multiple) site mutations. In particular, the SARS CoV-2 spike glycoprotein that exhibits interaction between its receptor-binding domain and the human angiotensin-converting enzyme 2, and has been seen to mutate rapidly with devastating consequences, was modeled and mutated in silico. The intramolecular HINT scores for the alpha, beta, gamma, delta, and omicron variants confirm that all mutated trimeric spike protein structures are similarly or more stable than the wild-type. In addition, these scores show that the receptor-binding domains of these mutants are more stable than the wild-type. The HINT intramolecular scoring is very rapid, and presuming the structures are well-modeled is suggestive that it may be useful as a pre-screening tool for as-yet undiscovered mutants for any protein.

The concept of interaction networks has been recognized for decades as it is a way to systematize protein structure, especially regarding hydrogen bonding. Our work with HINT has continually highlighted the parallel and often more important hydrophobic interactions. For example, it is not a coincidence that the α1β1 (and α2β2) interface of hemoglobin are largely unaffected by the deoxy to oxy hemoglobin transitions (and dominated by hydrophobic contacts between the two subunits), while the α1β2 interface, with just as many hydrogen bonds, has few hydrophobic interactions (Abraham et al., 1997) and shifts significantly in the transition. Also relevant are the “unfavorable” hydrophobic interactions (including those involving water), which are similarly key elements of protein structure. To systematize all of these observations and to create a platform for structure modeling, we developed a 3D interaction mapping paradigm wherein all interactions (polar and hydrophobic, favorable, and unfavorable) between a residue and its structural environment are visualized as contours with specific colors, strengths, and positions for each residue in a protein (Ahmed et al., 2015; Ahmed et al., 2019). To control backbone conformation, we adopted a chess square schema in which the Ramachandran plot is binned into 8 × 8, 45° × 45° chess squares (Ahmed et al., 2015). Further binning or parsing by sidechain conformations were applied as follows: χ ₁ for asparagine, aspartic acid, cysteine, histidine, isoleucine, leucine, methionine, phenylalanine, proline, serine, threonine, tryptophan, and tyrosine, and χ ₁ and χ ₂ for arginine, glutamine, glutamic acid, and lysine. The former, thus, have three parses per chess square and the latter nine (AL Mughram et al., 2021). These maps can, interestingly, after binning, be clustered into limited sets of interaction motifs characteristic of residue type, solvent accessibility, and structural role. The number of clusters obtained is dependent on the population of the associated bin and the complexity of observed and anticipated interactions. For example, alanine sidechains show about four unique maps per chess square, while arginine can show as many as 18 maps in some of the nine χ ₁/χ ₂ parses. Altogether, there are about 18,000 maps for our soluble proteins data set—each represents a unique constellation of interacting atoms for that residue type, backbone conformation, and sidechain parse. For illustration, Figure 7 sets out example contoured maps for six diverse residue types, all in the α-helix region of the Ramachandran plot. Importantly, it is not the identity of those interacting atoms but their specific properties, e.g., hydrophobic, polar, hydrogen bond donor or acceptor, etc., and water molecules can fill these roles just as other amino acid residues, cofactors or small molecules can. Because these roles are agnostic in terms of sequence, we call this paradigm “3D interaction homology” to emphasize that it is not the residue identity but its interactions that drive the structure.

With this extensive set of in-hand data, ∼18,000 maps abstracted from ∼750,000 residues in our dataset, we have a complete, quantifiable picture of the hydropathic valence of each residue type and its potential roles in a protein’s hydropathic interaction network. Exploiting this set of maps involves “matching” each backbone-aligned map’s encoded interactions between residues. Because there are limited sets of such maps for residue type and backbone conformation, this is a comparatively simpler problem than de novo atomistic minimization. More importantly, however, using HINT scoring and its core interaction model has enabled our view of the structure to be built upon a free energy framework and to account for some more subtle features of interactions that are not necessarily detected by molecular mechanics-based approaches. For example, sidechain maps of aromatic residues, phenylalanine, tyrosine, and tryptophan, show evidence of pi–pi stacking and pi–cation interactions (AL Mughram et al., 2021). Also, the computational titration algorithm was applied to aspartic acid, glutamic acid, and histidine map generation and allowed pH tuning of those residues’ interaction profiles (Herrington and Kellogg, 2021). Likewise, an understanding of the large differences in structural role between the serine and cysteine was explored with our HINT/interaction maps approach (Catalano et al., 2021). While our first goal was to document the interactions in which these residues engage, we were also able to develop an improved understanding of their roles in protein structure: 1) serine is considered to be somewhat more polar than cysteine but is significantly more solvent exposed; 2) while both have similarly consistent interaction roles (∼50–60% favorable polar, Figure 8) regardless of their accessibility, clearly very few cysteines are found on the outside of proteins; 3) it is also interesting (Figure 8) that bridging (-S–S- bonded) cysteines exist in much more hydrophobic environments than their unbridged analogs, which may be mechanistically suggestive. In another more recent study, the aliphatic hydrophobic residues were examined in both soluble and membrane proteins (AL Mughram et al., 2023), which surprisingly revealed somewhat modest differences in interaction profiles for these residues in soluble proteins, where they are often buried, and membrane proteins, where they are exposed to lipids. The latter is particularly important because very few X-ray crystal or even cryo-EM structures of membrane proteins retain their native lipids (Guo, 2020).

In the last few years, AlphaFold2 and RoseTTAFold (Baek et al., 2021; Jumper et al., 2021) achieved an astonishing power of prediction of protein three-dimensional structures thanks to the application of deep learning methods that revolutionized the approach applied until then (Goulet and Cambillau, 2022) This new capacity has changed the landscape of structural biology and drug discovery in surprisingly fundamental ways. There remain many proteins poorly predicted by these tools, and the emerging information concerning membrane proteins has likely not been effectively incorporated in these algorithms, as those structures themselves are not well-understood, especially with respect to the role(s) of their native lipids. In some cases, the value of experimental structural characterization has even been questioned! It should be noted, however, that even if computational predictions generate accurately-folded protein models, they are not likely to produce models that have sidechain conformations accurate enough for structure-based drug discovery and design (Tong et al., 2021). Indeed, it is the subtle features of the structure that drive toward success in these endeavors. HINT and other modeling tools that focus on the subtle interactions beyond gross folding and salt bridges will continue to be relevant and useful.

The analysis of energetic interactions involving DNA by using HINT was exploited for the first time in the study of the effects of antineoplastic drugs, such as chlorambucil (an alkylating agent), and anthracycline antibiotics, such as doxorubicin, daunorubicin, and others (which act as intercalating agents). In the first example, HINT analysis of a computational model of the adducts derived by administration of chlorambucil to the shuttle vector plasmid pZ189 demonstrated that the methylenes, as well as the phenyl ring of this drug, can form favorable hydrophobic interactions with nucleotides near the adduct site in the minor groove of DNA, promoting alkylation at the N-3 position of adenine (Wang et al., 1994).

Next, HINT was used to study the selectivity of doxorubicin intercalation and binding in all 64 unique base pair quartet combinations. The results showed that the interactions between doxorubicin and the base pairs above and below the intercalation site are mainly polar and favorable, resulting from acid–base interactions between the heteroatoms of the antineoplastic drug and the nucleotide bases. In addition, favorable hydrophobic contributions are also present. Specificity was mainly associated with hydrogen bonds, in particular with a base pair three positions away from the intercalating site (Kellogg et al., 1998). These analyses were subsequently extended to explore the binding of six different anthracycline antibiotics (doxorubicin, daunorubicin, hydroxydoxorubicin, 9-dehydroxydoxorubicin, adriamycinone, and daunomycinone) to 32 different DNA octamer sequences. The analysis showed that the differences in free energy among the various compounds were in line with that experimentally observed, indicating that HINT was able to pinpoint the energetic contributions of ligand functional groups most relevant to sequence specificity (Cashman et al., 2003). A further investigation of 65 doxorubicin analogs and their complexes with eight octamer DNA sequences allowed prediction of the net energetic contribution of several functional groups in the tested ligands and detection of selective ligands for the different combinations of DNA sequences simulated in this study (Cashman and Kellogg, 2004).

Analysis of the interactions of the two antibiotics gentamicin and paromomycin with 12 designed analogs with ribosomal RNA also showed that rings III and IV of these compounds are involved in important polar interactions with rRNA (Cashman et al., 2001). Therefore, HINT has proven to be a sensitive tool for dissecting interactions between nucleic acids and their ligands and for assessing the strength of interactions with existing drugs, as well as for predicting possible modifications of these molecules to improve their affinity and/or selectivity.

HINT was tested for its ability to evaluate interactions of much more complex systems, such as protein–DNA complexes. An initial study was conducted on the interactions between estrogen receptors (ER) alpha and beta and the corresponding estrogen-responsive elements (EREs) near estrogen-regulated genes (Marabotti et al., 2007). By analyzing the structure of the DNA-binding domain (DBD) of ERα and the homology-based model of ERβ bound to the ERE sequence, it was possible to identify the residues that contributed most to the binding affinity. Furthermore, by mutating each nucleotide pair in the two halves of the ERE binding site with all other possible pairs, it was possible to understand how mutations in the different positions of ERE could affect the binding affinity in both complexes. The results showed that, consistent with the experimental results, ERα binds the consensus ERE sequence with higher affinity than ERβ and that few amino acids and bases of the consensus sequences are involved in specific interactions. Specifically, HINT was able to discriminate with high sensitivity the affinity of ERα/ERβ DBDs for ERE sequences, as well as for non-ERE sequences used as negative controls (glucocorticoid- and progestinic responsive elements), whereas DDNA (Zhou et al., 2005), another predictor of protein–DNA interaction energies available at that date, did not. We hypothesized that the reasons for this failure was that the DDNA predictor was based on a knowledge-based statistical potential trained on a reference database not including protein–DNA complexes. However, it is significant that the HINT code was not derived specifically from the analysis of DNA structures; therefore, this fact confirmed us the general validity of the hydropathic approach (Marabotti et al., 2007).

For this particular set of complexes, the specificity of protein–DNA sequence binding did not appear to be much affected by water molecules. However, given the importance of the contribution of water in the thermodynamics of DNA–protein recognition, this investigation was extended in parallel to 39 additional DNA–protein complexes for which a three-dimensional structure was available (Spyrakis et al., 2007b). Thus, it could be shown that the inclusion of water molecules at the interface between protein and DNA (the so-called “bridging waters”) in the energetic contribution calculated by HINT improved the correlation between its score and the experimental free energy of association, with a lower standard error. The fraction of bridging waters in this set of experiments was only 3.5% of the water molecules detected in the 39 crystallographic complexes, consistent with the percentage of water mediating recognition between proteins and DNA identified previously (Reddy et al., 2001). It was also possible to observe that the orientation and binding strength of these water molecules depended more on the nature of the amino acid sidechain than on the type of DNA bases.

A more comprehensive study on energy-based prediction of the specific recognition between amino acid residues and nucleotide bases was subsequently performed on a dataset of 100 high-resolution protein–DNA complexes (Marabotti et al., 2008) (Figure 9) and used to predict specific contacts between amino acids and nucleotide bases in a set of 45 zinc finger–DNA complexes identified by phage display selection (Ghosh et al., 2006). By applying the HINT code, three main questions could be answered: i) Which amino acid–base pairs are energetically most relevant to achieve a specific interaction? ii) Are there energetic propensities that justify specific recognition of a nucleotide base by an amino acid? iii) Are bridging waters able to influence the specificity of amino acid–base recognition? The results showed that the amino acids that interact most frequently with nucleotide bases are Arg, Asn, Lys, Gln, Thr, Ser, Asp, and Gly; in fact, HINT calculated that the sum of their contacts accounts for more than 70% of the total number of contacts. Arg-G, Asn-A, Asp-C, Gln-A, Glu-C, and Lys-G appear to be the most energetically favorable contacts (Arg-G being the interaction that accounts for about 2/5 of the total HINT score for the complexes), while Asn-T, Asp-G, Gln-T, Glu-G, Ile-T, Leu-T, Met-T, and Val-T are unfavorable interactions. The analyses also showed that the same amino acid–nucleotide base pairs relevant to protein–DNA interactions are also particularly involved in water-mediated interactions.

For some amino acid–base pairs, it was also possible to calculate a “water enhancement factor,” that is, the ability of bridging waters to enhance the energetics of the amino acid–base interaction (Figure 10). Finally, based on the HINT score extracted from this analysis, it was possible to correctly predict more than 70% of the experimentally observed amino acid–base pairs in the zinc finger–DNA used as a test set. This percentage increased to nearly 90% when a relevance-weighted success descriptor, considering the relative energy relevance of each amino acid–base pair to the total protein–DNA recognition energy, was included. In this way, it was possible to show that amino acid–nucleotide base preferences could be explained by the energy-based analysis performed by HINT better than through qualitative approaches based on purely geometric considerations. Moreover, HINT also made it possible to predict unfavorable interactions, some of which are surprisingly well-conserved and usually involve the methyl group of thymine. This finding is interesting in that one might speculate that the amino acid–base interaction evolved before DNA development and was later adapted to DNA, but the presence of the thymine methyl group still continues to be a disruptive element in protein–DNA interaction (Marabotti et al., 2008).