Enthalpy-Entropy Compensation in Biomolecular Recognition: A Computational Perspective.

This mini-review provides an overview of the enthalpy-entropy compensation phenomenon in the simulation of biomacromolecular recognition, with particular emphasis on ligand binding. We approach this complex phenomenon from the point of view of practical computational chemistry. Without providing a detailed description of the plethora of existing methodologies already reviewed in depth elsewhere, we present a series of examples to illustrate different approaches to interpret and predict compensation phenomena at an atomistic level, which is far from trivial to predict using canonical, classic textbook assumptions.

Introduction

Enthalpy–entropy compensation (H/S compensation) entails a linear correlation between enthalpy and entropy changes in chemical processes where closely related structures or conditions are involved (e.g., reactions involving molecules differing only by a few functional groups, performed in different solvents or catalyzed by variants of the same enzyme) and is a general phenomenon affecting several physicochemical processes. From the point of view of biomolecular chemistry, it is particularly relevant in the fields of molecular recognition and drug design.[1] The search for new drugs and therapies often requires matching a candidate molecule with its target in order to stabilize as much as possible the resulting complex; in other words, the aim is to maximize the binding interaction of the candidate drug with a biological receptor (e.g., a protein or a nucleic acid). The strength of this interaction is evaluated through the binding free energy ΔGb, which is the main parameter to optimize. The process of drug design is iterative and has a strong structural component: crystallographic structures of the target in complex with a known ligand allow identification of the binding site(s) as well as a static representation of the ligand–target binding interactions. Based on this information, an initial virtual screening on a large library of compounds is usually performed to identify compounds capable of forming a stable complex, i.e., with a large negative value of ΔGb, with the target. The binding free energy for a ligand–target complex at a given temperature may be written as ΔGb = ΔHb – TΔSb where ΔGb is the sum of an enthalpic (ΔHb) and an entropic (TΔSb) term. Comparing the binding free energies of a family of compounds that bind onto the same target, H/S compensation may occur when structural differences among related ligands affect ΔHb and TΔSb in the same direction and to a similar extent, in such a way that the net effect of these enthalpic and entropic terms on ΔGb is negligible[1] (Figure ).

Figure 1

(a) Binding enthalpy (ΔHb, blue bars), entropy (−TΔSb, yellow bars), and free energy (ΔGb, green bars and black curve) measured experimentally by ITC for a series of ligands targeting HMG-CoA reductase. Data taken from ref (1). Despite the large variations observed in both enthalpic and entropic terms, the binding affinity remains nearly constant across the whole data set due to H/S compensation. (b) Plot of enthalpy (ΔHb) vs entropy (TΔSb) terms showing the linear relationship between them.

As it will be demonstrated in the following sections, H/S compensation phenomena cannot be explained by models of general validity but rather on a case-by-case approach. Hence, computer simulation plays a fundamental role in providing the means for connecting macroscopic thermodynamic binding signatures with molecular structure and dynamics, not only to explain experimental observations but also to design efficient binders with biomedical applications. In this mini-review, we will focus on examples of biomolecular recognition through noncovalent binding, presenting a series of examples to illustrate how computational techniques allow interpretation of the experimentally determined thermodynamic signatures of the H/S compensation processes. The scope of this work is not to provide a comprehensive overview of the computational methods for the calculation of binding free energies, which have been already covered in detail, but rather to introduce the reader to the complexity of the H/S compensation event by combining experimental and computational evidence. Thus, emphasis is placed on the use of different computational tools for interpretation of observed phenomena rather than on rigorous methods for predicting accurate ΔGb values.

Computational Approaches to ΔGb

Methods commonly used for ΔGb estimation can be broadly classified according to their accuracy and computational cost (Figure ). They can be roughly classified into: (i) equilibrium methods, such as free energy perturbation (FEP), thermodynamic integration (TI), and Bennett acceptance ratio (BAR), where the idea is to compute binding free energies through the structural perturbations between closely related states sampled with molecular dynamics (MD) trajectories; (ii) nonequilibrium methods, where the binding partners are progressively separated by an applied potential (steered molecular dynamics, SMD) and Jarzynski’s equality is used to reconstruct the free energy profile from information on the pulling force; and (iii) end-point methods, such as molecular mechanics combined with the Poisson–Boltzmann or generalized Born surface area continuum solvation (MM/G(P)SA), and linear interaction energy. A fourth class is represented by docking methods. The first two classes of methods yield a direct estimation of ΔGb and are vastly more computationally demanding than the third and fourth ones, which estimate ΔGb as a sum of terms and are significantly faster and less accurate. One of the main factors limiting the accuracy of these faster methods is the evaluation of the entropic contribution ΔSb, which is not straightforward. Given the wide use of these computationally efficient methods in different biochemical applications, we will summarize some of the approaches proposed to improve their estimation of the entropic contribution to binding.

Figure 2

Experimental and computational techniques commonly used in molecular recognition studies. Crystallographic structures are at the interphase of both sets of approaches and enable computational studies. Among the plethora of available computational methods, MM/PBSA and docking allow explicit modeling of H/S components to binding.

Molecular docking is at the core of drug design and allows screening a large number of candidate drugs onto the binding site of a receptor. Its computational efficiency lies in the use of scoring functions, which provide a means of ranking candidate ligands with a simple additive scheme, although with limited success in computer-aided drug design and development. In a recent report, Winkler explores the entropy terms employing several docking suites,[2] which normally enter the score with blunt approximations: ligand conformational entropy, originating from a multiplicity of conformations accessible to the ligand, is usually accounted for as a function of the number of rotatable bonds, while the loss of ligand translational and rotational entropy upon binding is computed as a function of the ligand molecular weight. From this study, it is apparent that even crude approximations to the binding entropy can improve the accuracy of ligand ranking.[2]

The popular MM/G(P)BSA method computes binding energies from classical molecular dynamics (MD) simulations using an additive scheme. ΔGb is estimated using a collection of structures extracted from the MD trajectory: explicit solvent molecules are removed, and binding enthalpy is calculated as the sum of an “internal” energy—an association energy based on the force field terms—and a solvation term calculated with an implicit solvent model. Values are averaged over the snapshots extracted from the molecular dynamics simulations, and in a first approximation the entropic contribution to binding is not accounted for. The success of these models lies in their ability to capture at least partially the flexibility of the binding interface, which is harder to represent with docking methods. Inclusion of entropy effects is the bottleneck for this type of calculation. The variation of conformational entropy upon binding can be computed through the normal-mode analysis of the snapshots. When such a term is included in MM/G(P)BSA calculations, it essentially dominates the calculation cost; for this reason, more approximate and computationally efficient schemes have been devised.[3] A first option is to perform the normal-mode analysis on a truncated system; entropies calculated with this approach are a good approximation to the full normal mode calculation, but care must be taken in the choice of dielectric constant used for optimization and frequency calculation, as it affects ΔSb. The interaction entropy approach is an alternative method that estimates ΔSb directly from the MD trajectory without further calculations through an exponential average. This method is less expensive and at least as accurate as normal-mode analysis on truncated systems but may suffer from errors deriving from the numerical instability of the exponential average, more prone to error propagation. To cure this problem, a truncated cumulant expansion can be used.[4] Quasi-harmonic analysis of MD trajectories, which is based on computing the covariance matrix of atomic coordinates, has shown serious convergence problems. A different approach involves calculation of the conformational entropy term from solvent-accessible and buried surface areas.[5] A further contribution to the binding entropy results from the restriction of the external degrees of freedom (translation and rotation) of the ligand upon binding to the receptor. Correction schemes based on the accessible volumes have been proposed to account for this effect.[4]

Experimental Measurement of H/S

Experimental ΔGb associated with biomolecular interactions (Figure ) is usually determined by isothermal titration calorimetry (ITC). With this technique, ΔHb and ΔGb are determined independently from heat measurement, and ΔSb is then obtained by subtraction. While ITC provides much more robust and reliable results than the older van’t Hoff analyses, which rely on the measurement of dissociation constants within a temperature range, still ΔHb and ΔSb values are not independent of one another, which may lead to artifacts. Nuclear magnetic resonance (NMR) can help unravel the thermodynamic signature of macromolecular binding phenomena complementing affinity data with structural information with a series of techniques, including transferred nuclear Overhauser effect (trNOE), saturation-transfer difference (STD), chemical shift perturbation (CSP), and relaxation experiments. Bio-layer interferometry (BLI) offers an alternative strategy to analyze biomolecular interactions: it is based on the measurement of interference patterns from white light reflected by two surfaces, an immobilized ligand and an analyte (in solution). This technique provides information not only on binding affinities but also on kinetic rate constants and allows quantitation of the analyte. A different class of biosensors is based on surface plasmon resonance (SPR). These biosensors exploit the sensitivity of a plasmonic material to the refractive index of its surroundings. A specific molecular receptor is immobilized on the material surface, allowing binding of the analyte. Binding involves a change in the refractive index on the sensor surface which elicits a signal in the plasmonic material. Several studies have applied this technique to investigate molecular interactions including determination of binding affinities, enthalpies, and entropies.

H/S Compensation and Interaction Strength

Central to biomolecular recognition, the extent of H/S compensation has been shown to be a function of the interaction tightness.[6] For extremely weak interactions, such as van der Waals complexes, ΔHb is small and varies slowly with structural modifications, while a strong entropic penalty arises from the loss of external degrees of freedom. Intermediate interaction tightness applies to the majority of ligand binding and protein–protein interaction events mediated by different types of hydrogen bonds, salt bridges, and van der Waals interactions, ΔHb ≈ TΔSb; it is in this range of binding affinities that H/S compensation can be observed as the two terms have approximately the same weight. For extremely tight binding, as can be the case for covalently bound drugs, ΔHb > TΔSb, and H/S compensation is no longer observed.[6]

On the Physical Origin of H/S Compensation

Several explanations have been proposed for this phenomenon, whose very existence has been questioned, and often attributed to an artifact arising from errors in the experimental determination of ΔHb and ΔSb values. Possible explanations invoke solvent structure, hidden Carnot cycles or a consequence of finite specific heat capacity, multiple weak interactions, quantum confinement, and limited free energy windows. Interestingly, it has been suggested that H/S compensation may be of evolutionary and functional advantage, providing a thermodynamic homeostasis that prevents harsh changes in free energy profiles. The intense debate on the nature of the H/S compensation phenomenon in the theoretical and computational community resulted in a rich literature that attempts to unravel this complex problem with simplified models.[7] A seminal theoretical work by Ryde addresses a fundamental question in H/S compensation regarding noncovalent binding, i.e., whether contacts dominated by hydrogen bonds or van der Waals forces can present different H/S compensation profiles owing to the fundamentally different nature of the weak interaction involved.[8] The report discusses this by using model systems fully dominated by either interaction type and considering the effects of molecular size and solvation. The general conclusions are that no fundamental difference is observed for different noncovalent interaction types in the H/S balance, leading to compensation and obtaining comparable ΔHb vs TΔSb linear relations. Of note, this linearity becomes blurred when solvent molecules are added to the models. This analysis determines the two main players affecting H/S compensation in molecular recognition: (i) the flexibility and roto-translational freedom of the binding partners, which are deeply affected by complex formation, and (ii) solvent effects. In the following sections, we will provide real-case examples in which both sets of factors can dominate binding affinity in a system-dependent manner.

Robustness vs Flexibility

Conformational Entropy Terms

Conformational or configurational entropy is classically associated with the reduced number of translational, rotational, and vibrational degrees of freedom of the ligand upon binding and contributes unfavorably to the binding energy. Coupling between these degrees of freedom makes it difficult to split this contribution into additive terms.[9] Historically, acknowledgment of the importance of flexibility effects has fundamentally changed the interpretation of macromolecular recognition phenomena. The “lock and key” model was first proposed by Fischer in 1899 as a way to interpret enzymatic activity; this model bases recognition on structural complementarity, i.e., the ligand and receptor having predefined conformations that fit with each other, in such a way that electrostatic and van der Waals interactions are established between the binding partners at a fixed geometry. This view, which implies representation of the binding partners as rigid entities, implicitly emphasizes the role of ΔHb at the expense of ΔSb. Over the last years, this model has been replaced by the so-called “hand and glove” model, which in turn assumes a certain flexibility of the binding partners and a mutual conformational influence, accounting for both enthalpic and entropic contributions to binding.

Despite being just one component of ΔSb, the ligand’s conformational entropy can affect binding affinity even when small structural modifications are involved. A very clear example of this effect was recently reported by Ernst and co-workers.[10] By designing ligands for the mannose-specific bacterial lectin FimH, they compare the binding affinities of the natural ligand, n-heptyl α-d-mannoside, and a seven-membered ring mannose mimetic. The two molecules only differ by a single bond in the ring. The larger analogue shows an approximately 10-fold lower affinity than the natural ligand. ITC experiments reveal that this difference in affinity essentially originates from the entropic term. Further ITC experiments conducted at higher temperatures allow a partial decoupling of the solvation and conformational terms of this entropic difference. This is crucial because it is extremely difficult to separate conformational and solvation contributions from ITC data. In this case, an essentially equal solvation entropic contribution is observed for the two ligands, which is to be expected given the similar structure and distribution of polar groups of the two compounds. The high structure similarity of the two ligands also implies in this case a similar entropic contribution from the lectin receptor upon binding, as revealed by NMR CSP experiments. The molecular origin of this penalty is revealed by metadynamics MD of the free ligands: the additional bond in the seven-membered mimic imparts a much larger conformational flexibility than that of the mannopyranoside, with multiple isoenergetic shallow minima. Insertion of the ligands in the binding pocket involves a strong rigidification, meaning that little residual flexibility is maintained in the binding pocket, with a larger conformational penalty for the more flexible ligand. The importance of this elegant work lies in its decomposition of the entropic contributions to binding in its constituting terms, which allows the molecular interpretation of the thermodynamic binding signatures and directs the computational study to the detailed analysis of the conformational profiles of the unbound ligands. This example highlights the importance of the residual flexibility of the bound ligand: while the most intuitive strategy in a drug design problem is to maximize the number of ligand–receptor interactions, improving the enthalpic contribution to the binding energy—and worsening the entropy term—less attention is usually paid to improving the entropy term. A similar impact of ligand rigidity on binding affinity was recently reported by Jiménez-Barbero and co-workers in the recognition of histo-blood group antigens A and B vs the more flexible N-acetyl-d-lactosamine (LacNAc) by human galectin-3.[11]

The paradigm shift in the interpretation of biomolecular recognition (“lock and key” versus “hand and glove” models) has also led the community to enquire into the relationship between the robustness of noncovalent interactions and binding energy. Barril and co-workers investigated a large set of hydrogen-bonded protein–ligand complexes using SMD simulations to calculate the work necessary to break individual interactions.[12] These values are a measure of the robustness of the given interaction. Results reveal that a high content of robust hydrogen bonds is rare and that most complexes feature a single robust contact surrounded by a looser network of other interactions. These results indicate that successful ligand–receptor pairings respond to H/S compensation by attaining a flexible binding mode that optimizes the trade-off between enthalpic (“robustness”) and entropic (“residual flexibility”) contributions. In some cases, no robust interaction is present, questioning the validity of drug design approaches based on individual, well-defined binding modes.[12] Taking this idea to the extreme, Borgia et al.[13] reported an ultrahigh affinity protein–protein complex which completely bypasses the requirement of structural complementarity and a defined binding mode; despite showing a picomolar affinity, the complex fully retains the disorder of the two binding partners (i.e., undergoes no entropic penalty).

A special instance of “robustness” versus “residual flexibility” in ligand binding is the cooperativity of noncovalent interactions and its implications for H/S compensation. Such cooperativity can be negative (i.e., the combined interaction yields a smaller binding free energy than the sum of the individual binding free energies) or positive (when the combined ΔGb is greater than the sum of the individual terms). Two main models have been proposed to interpret cooperativity phenomena: partially bound states[14] and structure tightening.[15]

The partially bound state model, proposed by Hunter and co-workers,[14] postulates that in systems regulated by multipoint interactions the bound state can be represented as an ensemble of different complexes, each with a different distribution of contacts among all possible ones (Figure a). Of all possible states, one is fully bound (i.e., all possible contacts are formed), while in the others the partial loss of contacts is compensated by a larger flexibility resulting in a dominant entropic benefit. Although conventional, fixed-charge force-field MD simulations are suited to analyze multivalent binding behavior as they allow sampling the conformational space of bound complexes, cooperativity is elusive to these methods due to their inability to account for the polarization of interacting groups.

Figure 3

(a) Partially bound state model for multipoint interactions. The increasing number of possible contacts in larger ligands increases the multiplicity of partially bound states; favorable ΔΔSb exceeds ΔΔHb penalty from reduced contact frequency. (b) Structural tightening model for multipoint interactions. The increasing number of contacts shortens noncovalent interactions (positive cooperativity), and favorable ΔΔHb exceeds ΔΔSb penalty.

According to the structural tightening model, proposed by Williams and co-workers,[15] positive cooperativity observed for multipoint interactions has an enthalpic origin: noncovalent bond distances are shorter in multipoint interaction complexes than in complexes where the same interactions appear individually (Figure b). This originates from an increased organization of the bound state when multipoint interactions are established, which translates into in a favorable enthalpic stabilization that exceeds the entropic penalty. Quantum mechanical calculations are capable of reproducing the mutual polarization of functional groups involved in cooperative interactions, although commonly from a static perspective. An example of computationally studied positive cooperativity that fits into the structural tightening model is the binding of biotin to avidin and streptavidin.[16] When interacting with both receptors, biotin forms three hydrogen bonds. The contribution of these hydrogen bonds to the binding affinity is analyzed on cluster models with a combination of density functional theory (DFT) and Møller–Plesset quantum mechanical methods, uncovering a significant positive cooperativity. Calculations unveil the importance of a key aspartate residue, which mediates the cooperativity effect; this aspartate is directly in contact with biotin through a hydrogen bond in streptavidin and is one residue away in avidin. This results in an ∼4 kcal mol–1 larger cooperativity effect to binding in streptavidin than in avidin. Modeling shows that the complete removal of the aspartate from the avidin binding site further disrupts cooperativity with marked increase of the remaining hydrogen bond lengths.

Enzymatic Activity Benefits from Ligand-Induced Conformational Entropy Variations in the Receptor

Until now, we have discussed the importance of ligand entropy contribution to binding free energy, but also receptor entropy can tune molecular recognition and regulate cooperativity. One example that stresses the critical role that can be played by a receptor’s flexibility was recently published by Veglia and co-workers.[17] They used a combination of ITC, NMR spectroscopy, and MD simulations to investigate the allosteric cooperativity in the catalytic cycle of cAMP-dependent protein kinase A. Results show that cofactor (ATP) binding rigidifies the protein, contributing an unfavorable entropic term. This initial rigidification facilitates substrate binding, which in turn increases enzyme conformational dynamics (i.e., favorable entropy contribution) and provides an overall positive cooperativity. Conversely, negative binding cooperativity is observed for the phosphorylated product and ADP, which facilitates the rate-limiting ADP release step. MD simulations allow decomposition of the configurational entropy term on a per-residue basis. Mapping the per-residue configurational entropy on the enzyme structure reveals that entropic changes are not localized but rather distributed over the enzyme structure. Overall, these results show that the protein’s conformational entropy is involved at all stages of the catalytic cycle and can fine tune ligand affinity.

The balance between enthalpy and entropy affects not only substrate affinity in enzymes but also their catalytic activity defined by the rate-limiting activation free energy. According to Eyring’s equation:The reaction rate constant kr depends exponentially on the activation free energy ΔG‡. The stark acceleration of enzymatic reactions compared to uncatalyzed ones results from deep changes in the thermodynamic signature of activation parameters. Binding into an enzyme active site primes the substrate for catalysis; moving the substrate from the solution bulk to the enzyme active site entails a loss of rotational and translational entropy which is “prepaid” by the formation of the Michaelis–Menten complex, providing an entropic advantage. Additionally, enzyme-active sites are highly preorganized to minimize the enthalpic component of the activation barrier by preferentially stabilizing (i.e., “binding”) the transition state according to Pauling’s paradigm of biocatalysis. Åqvist and co-workers proposed a scheme to calculate thermodynamic activation parameters of chemical and enzymatic reactions from MD simulations,[18] which not only affords ΔH‡ and ΔS‡ values in excellent agreement with experimental data but also provides valuable insights into reaction mechanisms. As an instructive example, we have selected an application of this computational scheme to cold-adapted enzymes.[18] These systems are particularly attractive due to the unusual partitioning of their activation parameters: compared to their mesophilic orthologs, cold-adapted (psychrophilic) enzymes show lower enthalpies and more negative (higher absolute value) entropies of activation. This evolved feature is also a H/S compensation phenomenon: lowering the activation enthalpy ΔH‡ decreases the activation barrier ΔG‡ but is compensated by the larger −TΔS‡ term. This larger entropic component protects the activation free energy from its natural damping with the decreasing temperature, allowing psychrophilic enzymes to achieve high reaction rates in cold conditions. This technique not only reproduces well the experimental activation parameters but also provides insight into the enthalpy–entropy balance of the two types of enzymes. Contrary to what may be expected, structural differences are not observed in the active site—which can heavily influence the reaction rate by interacting directly with the substrate—but are distributed on the protein surface. Overall, psychrophilic enzymes have mechanically softer protein–water surface interfaces than their mesophilic counterparts, which explains their increased entropic contribution to the free energy barrier.

H/S Compensation Originating from Drug-Resistant Mutations

Another example of how the receptor structure can change the thermodynamic landscape of binding was proposed by Schiffer and co-workers.[19] Mutations can alter the thermodynamic signature of binding, leading to a decrease in binding affinity, and represent a major challenge for drug design methods. Even if changes in the binding affinity are small, they can hide major changes between enthalpic and entropic terms, as shown in this example. Schiffer and co-workers focus on HIV-1 protease, a paradigmatic receptor for drug design studies, and consider how mutations affect the binding affinity toward a set of structurally diverse inhibitors.[19] The HIV-1 protease variant under study, called Flap+, shows a 1–3 kcal/mol reduction of the binding affinity for all six considered inhibitors compared with the wild type as measured by ITC; this comparably small difference between the two receptors is the consequence of a large H/S compensation phenomenon, revealing much larger (5–15 kcal/mol) variations in the ΔHb and TΔSb terms. MD simulations attribute this difference to a marked change in flexibility in specific regions of the protein, highlighting once more the importance of receptor flexibility and the dependence of H/S compensation phenomena on the specific features of the system under study.

Solvation

Water is key to biomolecular recognition and can exert its influence on association processes through a variety of different mechanisms.[20] In this section, we will discuss some examples of how solvation can determine the thermodynamic signature of association events, with focus on the insight that can be provided by computational methods.

Solvation is an Intrinsically Compensatory Process

As molecular recognition takes place in aqueous solution, modifications of the hydration shells of the binding partners upon binding also contribute to the thermodynamic signature of binding. However, while ligand–receptor interactions can be analyzed in terms of intuitive and computable concepts—shape and charge complementarity for enthalpy and flexibility for entropy—no such basic notion is available for the role of solvation. As an example, protein surface hydration is a complex and dynamic phenomenon resulting from a balance between bulk water and the propensity of exposed side chains to be solvated, which in turn is an intricate function of the protein structural features. When considering ligand–receptor binding from the solvent point of view, two fundamental processes take place: (i) desolvation of the ligand and the receptor’s binding cavity and (ii) solvation of the resulting complex in the binding region. Normally, water desorption is characterized by unfavorable enthalpy and favorable entropy—transfer to the solvent bulk increases the accessible number of degrees of freedom—while the opposite effects are associated with water adsorption (i.e., hydration). Additionally, polarization and hydrophobic entropies have been theoretically defined to account for changes in the orientational freedom of solvent dipoles upon ligand binding.[9] In both cases, entropy and enthalpy contributions go in opposite directions, which indicates that H/S compensation is an intrinsic feature of solvation and contributes to the total H/S compensation of binding. The tightness of water binding at a ligand–receptor complex surface and, consequently, the strength of its H/S compensation character depend on the rigidity of the complex itself, as highlighted by Crane–Robinson and co-workers.[21] Solvation amplifies the usual enthalpic gain and entropic loss associated with binding: rigid complexes present tighter first hydration shells than looser ones, adding to their enthalpic stabilization and entropic destabilization. Immobilization and release of water molecules is a ubiquitous phenomenon that always accompanies ligand–receptor recognition in solution. This phenomenon is now widely recognized but very difficult to tackle with computational methods owing to the large number of degrees of freedom involved in the treatment of solvent molecules. For this reason, solvation entropy is often neglected or treated with crude approximations that can compromise interpretation of experimental binding affinities.[21]

Hydrophobic Effect

The hydrophobic effect has long been the paradigmatic explanation for the association of mostly nonpolar binding partners interacting only through weak interactions. In this kind of complexes, a net release of water molecules upon binding provides a net entropic gain that favors association. In the classical view of this phenomenon, water release entails an enthalpic penalty as water molecules move from the highly structured receptor surface to the less ordered solvent bulk; this penalty, however, is thought to be overcome by the entropic gain, making water release a favorable process. In some cases, however, it has been found that displacement of solvent molecules can also be enthalpically beneficial. For instance, Cockroft and co-workers analyzed hydrophobic molecular association by decoupling solvent effects from van der Waals interactions, the latter being evaluated by ab initio and DFT quantum mechanical calculations and gas-phase measurements.[22] Discrepancies observed between interaction energies obtained in the presence and absence of water suggest that enthalpic, cohesive solvent–solvent interactions can be the major driving force for the association of nonpolar species in solution. Capturing such desolvation effects often requires explicit modeling of water molecules because a continuum solvent representation does not account for structured or semistructured water at the ligand and receptor surfaces.

Frustrated Solvent Local Structure

Depending on the topology of protein surfaces, adsorbed water molecules interacting with both the protein side chains and other waters can exhibit different entropic and enthalpic signatures compared to bulk ones.[23] In general, they show more favorable enthalpy and less favorable entropy than bulk water; this, however, is not always the case. Kurtzman and co-workers[23] analyzed the binding site hydration of six structurally diverse proteins using hydration site analysis and measures of local water structure through MD simulations. Their results showed that certain protein structures can adsorb water by providing a lower enthalpic stabilization to these water molecules compared to the bulk. Water molecules at such sites are thus “frustrated”, and their transfer to the bulk solvent upon ligand binding contributes favorably to both binding entropy and enthalpy (Figure ).

Figure 4

Energetics of water displacement upon ligand binding. While the entropic component is favorable and the enthalpic component is normally unfavorable, when frustrated water molecules reside at the receptor’s binding site, the process becomes also enthalpically favorable upon transfer to the bulk solution due to stabilizing solvent–solvent interactions.

Nontrivial Role of Charged Groups

Barril and co-workers used SMD simulations to show that hydrogen bonds involving charged groups are only slightly more robust (∼1 kcal/mol difference in mean work to break the interaction) than neutral hydrogen bonds.[12] This can be explained as a compensation effect: charged groups benefit from an additional electrostatic contribution, which strengthens the interaction, but the desolvation penalty of charged groups is higher than that of neutral ones.[12] In biomolecules, the formation of salt bridges between charged groups contributes to binding free energies with a distinctive signature. Indeed, solvent-exposed charged groups interact nonspecifically with environmental counterions to achieve global neutrality. Upon formation of a salt bridge between the ligand and the receptor, counterions are released to the bulk solvent in an enthalpically balanced and entropically favored process.[21] This net entropy gain suggests that increasing the number of ligand–receptor salt bridges can be a useful strategy to increase binding affinity

Halogen versus Hydrogen Bonding

Halogen bonding can be described as a highly directional net attractive intermolecular interaction between the electrophilic region (σ-hole) of a halogen and a nucleophile. Halogen bonds can be seen as analogues of hydrogen bonds, with a less polar character. One consequence of this larger hydrophobicity is a reduction of the solvation penalty upon binding. Ho and co-workers[24] systematically analyzed this effect and its repercussion on H/S compensation determining the crystal structures of DNA Holliday junctions in complex with halogenated uracil bases. In this study, they analyzed the different thermodynamic contributions to the binding affinity stemming from halogen bonds compared with the classic hydrogen bond interaction. Calculation of solvent-accessible surfaces is used to demonstrate that burying a halogen instead of a polarized hydrogen is favorable due to solvation effects. This provides a way to break H/S compensation: comparing the binding energies for a Br bond versus the corresponding hydrogen-bonded construct, a net gain in binding affinity is observed with both enthalpic and entropic stabilization.

Key Role of Desolvation in Saccharide Binding

Carbohydrates are relatively rigid ligands rich in hydroxyl groups which are often seen as structured water molecules; their rigidity is usually matched by their receptors’ binding sites, which tend to be highly preorganized. Furthermore, the abundance of polar groups entails a high desolvation penalty to their binding. The combination of these two effects leads to a significant coincidence in the positions of structural waters in the receptors’ binding site and those of the carbohydrate hydroxyl groups in the bound state, so that binding is enthalpically balanced. In this way, water displacement by carbohydrate hydroxyl groups provides a favorable entropic drive for binding. Ernst and co-workers[25] recently combined ITC, X-ray crystallography, and MD simulations to determine the binding affinity of six oligosaccharides that act as antagonists for PapG-II and shed light on the factors affecting the stability of the interaction and the role of water. These oligosaccharides show H/S compensation but with an unusual dependence on system size: smaller oligosaccharides show higher enthalpic gains and entropic penalties than larger ones. This result is opposite to the common size dependency trend in that larger ligands are expected to establish more enthalpically stabilizing contacts and undergo more significant conformational entropy restrictions than smaller ligands. However, MD simulations revealed that changes in binding energy along the ligand series cannot be traced back to the number of intermolecular interactions and restriction of the molecular flexibility alone, as all the ligands show very similar flexibility profiles. In this case, solvation effects dominate both enthalpy and entropy and become the main thermodynamic driving force: on one hand, the enthalpic cost associated with hydroxyl group desolvation overpowers the enthalpic gain of establishing ligand–receptor contacts, leading to a net enthalpic penalty that increases with the increasing ligand size and number of hydroxyl groups. On the other hand, the entropic term follows an opposite trend, becoming more positive (favorable) with the system size as a larger number of water molecules around the hydroxyl groups are released into bulk water upon binding.

Woods and co-workers reported careful experimental and computational efforts to increase the binding affinity of a trimannoside toward concanavalin A (Con A) by structural water displacement.[26] In a first work, the authors analyze through a combination of ITC, NMR, and MD simulations the interaction of Con A with two trimannoside ligands (one natural and one artificial) differing by a C2H4 unit. Where the natural ligand presents an OH at the C2 position, the artificial one presents a hydroxyethyl moiety. This structure difference responds to a precise design: a highly ordered water is present in the binding site that is not displaced by the natural ligand, but the additional chain introduced in the artificial one is expected to displace this water molecule, resulting in an improved binding affinity. ITC measurements report a lower binding affinity for the artificial trimannoside over the natural one but a more favorable entropic contribution to binding. This more favorable entropic term is initially attributed to structural water displacement, apparently corroborating the design hypothesis. A subsequent work by the same authors disproves it: the crystal structure of the complex of Con A with the artificial ligand reveals that the added chain does not displace the structural water molecule, even if it distorts its surrounding compared with the complex of the natural ligand. This structural information provides the basis for an accurate computational analysis of ligand binding through thermodynamic integration (TI). In this work, the reduced entropic penalty measured for the artificial trimannoside is attributed to an H/S compensation phenomenon, i.e., counterbalancing a reduced enthalpic contribution to binding compared with the natural ligand. A further analysis of this system relies on computation to assess the properties of the water to be displaced. The authors hence focus on the binding free energy of this structural water molecule to assess the thermodynamic profile for its association to both the bound and unbound form of Con A.[26] Interestingly, they find that the model used to represent water (TIP3P, TIP4P, or TIP5P) has a large effect on the lability of structural water molecules. This finding is extremely significant as a caveat for drug design workflows that include MD simulations including explicit solvent, as results may largely depend on technical aspects of the simulations.

Conclusions

The practical cases reviewed here show that, despite being a general phenomenon affecting a wide range of biomolecular interactions, H/S compensation has no obvious dependence on the system features. The dual nature of the entropic change that accompanies association processes (desolvation and conformational entropies) prevents drawing general conclusions on the effect of properties even as simple as ligand size since solvation and flexibility can alternatively gain the upper hand depending on the specific characteristic of each system. In particular, the examples presented in this work, which combine both experimental and computational observations, caution against beforehand assumptions on the role of entropy in determining the thermodynamics of binding processes, demonstrating that the entropy term, particularly when binding takes place in water, should never be neglected.

As the levels of complexity and resolution that can be tackled with simulations will keep increasing in the upcoming years, experimental and computational chemists must work hand-in-hand to base method development on firm experimental grounds. Only through collaborative, cross-discipline approaches, the accurate thermodynamics of underlying molecular recognition events, in which H/S is central, will be reliably and consistently unraveled. When achieved, this will dramatically increase the prediction capacities of computational methods to accelerate the discovery of more efficient medicines and therapies.