How Can Hydrophobic Association Be Enthalpy Driven?

Hydrophobic association is often recognized as being driven by favorable entropic contributions. Here, using explicit solvent molecular dynamics simulations we investigate binding in a model hydrophobic receptor-ligand system which appears, instead, to be driven by enthalpy and opposed by entropy. We use the temperature dependence of the potential of mean force to analyze the thermodynamic contributions along the association coordinate. Relating such contributions to the ongoing changes in system hydration allows us to demonstrate that the overall binding thermodynamics is determined by the expulsion of disorganized water from the receptor cavity. Our model study sheds light on the solvent-induced driving forces for receptor-ligand association of general, transferable relevance for biological systems with poorly hydrated binding sites.

Introduction

It is becoming widely recognized that water−the common environment for most biological processes−plays a significant role in binding (thermo)dynamics, being far more important than just a passive, embedding medium. This is particularly relevant for hydrophobic interactions, in which water-related effects are typically regarded as the origin of the entropic driving force. Nevertheless, our knowledge about how water structure and dynamics evolve during such binding events and what are the underlying thermodynamic contributions is still sparse.

One of the recent, key advances is the distinction among different length scales involved.[1−4] It is now generally accepted that in the limit of small solutes, e.g., small hydrocarbons, the surrounding water structure and hydrogen bonding are not considerably affected. In this case, hydration thermodynamics appears to be dominated by the entropic effect of restricting spontaneous solvent fluctuations to only those permitting the presence of the solute.[5−8] The resulting solvent-mediated interactions, exemplified by the well-investigated case of methane pair in water,[9−12] display a characteristic, entropy-stabilized free energy minimum for the direct contact of two solutes. However, the spontaneous assembly of small isolated hydrophobic molecules is not observed,(13) owing to a comparatively larger configurational volume of the solvent-separated pairs.

Hydrophobic-driven association usually involves at least one interacting partner carrying an extended nonpolar region. It affects the hydrogen-bonding network of interfacial water molecules, inducing a variety of large-length scale hydrophobic effects.[14,15] The actual solvent behavior critically depends on the strength of solute−solvent attraction and varies from the persistent hydration of strongly interacting objects,[16,17] through the formation of a thin vapor-like interface next to flat or moderately concave hydrophobic surfaces,[18−20] to the complete dewetting of sterically hydratable regions.[21,22] The importance of such effects for the assembly of nanoscopic bodies was investigated in a number of molecular dynamics (MD) simulation studies.[23−34] Attempts to characterize the underlying thermodynamic signatures in this context were, however, limited to relatively simple solutes, like plates or ellipsoids,[23,28,33] and support the conventional view of entropy-driven hydrophobic association.

Is this conventional picture directly transferable to biologically relevant binding scenarios? They typically involve a receptor cavity, whose degree of hydration may vary significantly depending on geometry and composition,(15) and a ligand, whose hydration belongs to the aforementioned small-length scale regime. Following an elegant reasoning by Dunitz,(35) it is usually assumed that the release of receptor-bound water is entropically favorable. This interpretation may hold even in the case of hydrophobic cavities, in which water is weakly bound. In such case, MD simulations indicated the formation of hydrogen-bonded water clusters,[21,36−38] whose displacement upon binding should lead to the expected entropy dominated hydrophobic association. On the other hand, fundamental insight from nuclear magnetic resonance experiments suggests that water molecules permeating protein cavities may be also disordered[39,40] with their displacement accompanied instead by an entropy loss,(40) in line with a recent MD study.(41) These important findings remain largely unexplored in the context of cavity−ligand association, possibly contributing to unexpected effects, like enthalpy-driven binding of nonpolar ligands to the highly hydrophobic pocket of mouse major urinary protein (MUP).[42−44]

Here, we investigate the driving forces for model nonpolar cavity−ligand association (Figure 1).[30,45] This type of hydrophobic binding is peculiar in that it simultaneously involves both small- and large-scale alternative hydration regimes. Using explicit solvent MD, we derive the potential of mean force (PMF) for cavity−ligand interaction at five different temperatures. From the temperature dependence of the free energy, we obtain a complete picture of thermodynamic signature profiles along the association coordinate, thus far limited to simple solutes. As both associating partners are devoid of internal degrees of freedom, we directly capture water thermodynamic contributions and relate them to the ongoing changes in solute hydration. To our knowledge this is the first study elucidating the interplay between enthalpic and entropic components in the context of model cavity−ligand binding. We propose an explanation on why the observed association appears enthalpy driven, contrarily to what is usually expected for different hydrophobic systems.

Figure 1

Snapshot and schematic representation of the explicitly solvated hemispherical cavity and spherical ligand (L) used in this study. Note that (ξ = 0) corresponds to the wall surface.

Methods

Molecular Model and MD Simulations

The hemispherical cavity of 0.8 nm radius was embedded in a rectangular paraffin-like wall (Figure 1), constructed as a hexagonal close packed (HCP) grid (lattice constant 0.125 nm) of pseudoatoms interacting through 6−12 Lennard-Jones (LJ) potential with parameters εp = 0.0024 kJ mol−1 and σp = 0.4152 nm. Note that despite a small εp value, the well depth for flat wall−water interaction potential is 2.93 kJ mol−1 due to tight packing of wall particles, and hence, the wall should be regarded as realistic hydrophobic material (see the Supporting Information for details). The ligand was modeled as one neutral LJ sphere (methane parameters: εm = 1.23 kJ mol−1 and σm = 0.373 nm(46)), and the TIP4P model was used for water molecules.(47) Pairwise LJ interactions were treated using standard mixing rules. The system used in MD simulations consisted of two identical 3.5 × 3.3 nm2 pocketed walls with adjacent ligands, mirrored along the Z axis such that they made two opposite sides of a 3 nm thick box filled with 1030 water molecules. Periodic boundary conditions were enforced with a box size equal to the system size in X,Y directions and 10 nm size in Z direction, i.e., with vacuum behind the walls. Association of each ligand with its pocket was sampled independently along ξ (the system symmetry axis) from 1.1 nm (in the bulk region) to −0.4 nm (inside the pocket) over 31 consecutive windows (0.05 nm apart) using the umbrella sampling method.(48) The corresponding PMF, W(ξ), was calculated using the weighted histogram analysis method.(49) In all simulations, the distance between the two ligands was at least 2 nm (i.e., when one of them was in the bulk region, the other occupied the pocket interior), and it was assured that their movements were not correlated. All simulations were carried out with the CHARMM software.(50) For additional details we refer to refs (30 and 45). Five sets of independent simulations were performed at T = 298, 308, 318, 328, and 338 K. Simulation for each umbrella potential window was preceded by a number of 1.1 ns equilibrating runs during which the wall separation was iteratively adjusted until water density at the center of the solvent box matched the experimental value at the given temperature and pressure of P = 1 bar(51) with tolerance ±2 g/L. The actual production runs (1 ns, preceded by 100 ps of equilibration) corresponded to NVT conditions with individual system volumes for each sampling window, in order to avoid barostat artifacts due to the presence of constrained walls (see ref (52) on this point). Owing to the system symmetry, effectively 2 ns sampling was obtained per ξ value, per temperature. Water density distribution maps were generated using xfarbe.(53)

Free Energy, Entropy, and Enthalpy Changes and Their Uncertainties

The Gibbs free energy for the system with the ligand at a given ξ value reads where G(∞) was set to 0 by shifting to 0 the average G for ξ ∈ [0.95, 1.1] nm. The corresponding entropy, S(ξ), can be obtained through the partial temperature derivative of the free energy. The corresponding enthalpy reads H(ξ) = G(ξ) + TS(ξ).

The water contribution to the total Gibbs free energy was obtained as GW(ξ) = G(ξ) − UCL(ξ), where UCL is the direct ligand−cavity interaction. Ligand−water, ULW, and cavity−water, UCW, interaction energies were calculated using force field terms with 1.2 nm cutoff as ensemble averages over MD trajectories at ξ values discretized into 0.01 nm bins. The water−water interaction energy was obtained as UWW(ξ) = H(ξ) − ULW(ξ) − UCW(ξ) − UCL(ξ).

The standard formula expresses—for a given system at constant pressure—the temperature dependence of G as a function of the system heat capacity, Cp. Assuming constant Cp for the considered temperature range, such equation can be rewritten in a parametric form as Differentiating with respect to T gives The results presented throughout this article correspond to S(ξ) numerical estimate for T = 298 K, obtained from optimal A and B values to fit eq 3 using five G(ξ) values.

Uncertainties on G and S estimates were determined as where X̃ is a value obtained using all data points, while X corresponds to N = 5 independent simulation blocks. Uncertainties δX for G, H, and −TS profiles along ξ are reported as vertical error bars (Figures 3 and 4).

Figure 3

(Thermo)dynamics of pocket hydration. (A) Average cavity occupancy along ξ. is the average number of water oxygens located at Z < 0. (B) Free energy and its entropy component as functions of pocket occupancy, N (note different energy scale for each profile).

Figure 4

Thermodynamic contributions along the binding coordinate. (A) Relative Gibbs free energy, G (red), enthalpy, H (blue), and entropic term, −TS (green), and their uncertainties (vertical bars; eq 5). (B) Water contribution to the Gibbs free energy, GW (orange), and decomposed interaction energies: cavity−ligand, UCL (thin black), ligand−water, ULW (green), cavity−water, UCW (black), and water−water, UWW (cyan).

We note that the uncertainties for entropic and enthalpic components for T = 298 K are larger than for T = 318 K, the midpoint of the considered temperature range (see the Supporting Information). Nevertheless, as no qualitative difference in the resulting signatures is observed at the two temperatures, we focused on the thermodynamic data for T = 298 K in order to remain in line with the majority of other studies conducted at this temperature. We stress that all conclusions of this study regarding the cavity−ligand binding process are independent of this choice.

Results

Hydrophobic Hydration

First, let us focus on the hydration of the ligand and the binding pocket when they are far from each other. Two clearly distinguishable hydration shells are formed around the ligand, in agreement with the expected picture for small-length scale hydration, as inferred from the water density map for ξ = 1.1 nm (Figure 2). In contrast, average water density inside the pocket gradually vanishes, indicating that dewetting takes place (Figure 3A). This can be attributed exclusively to the concave topography of pocket walls, as the water potential energy due to interaction with the solutes is lower at the pocket bottom (−4.6 kJ mol−1 minimum) than around the ligand (−0.9 kJ mol−1 minimum).

Figure 2

Water density distribution maps for key snapshots along ξ. Color coding is normalized such that ρ* = 1 corresponds to bulk water density of 998 g/L.

A previous analysis of pocket hydration revealed that the observed average water density results from intermittent expansions and retractions of the liquid phase rather than uncorrelated diffusion of individual water molecules characteristic of a “vapor-like” phase.(52) In order to quantify the thermodynamic effects due to pocket hydration, we calculated the system free energy, F, as a function of pocket occupancy, N, using the relation F(N) = −kBT ln P(N), where N is the number of water molecules with an oxygen atom located at z < 0, P(N) corresponds to its probability distribution obtained from MD trajectories with the ligand at ξ ≥ 10 nm, T is the temperature, and kB is the Boltzmann constant. The F profile reveals series of metastable states extending from N = 0 (empty pocket) through N ≃ 6 (average occupancy) to N = 11 (close to bulk-like density), indicating no apparent free energy change for pocket dehydration (Figure 3B). The corresponding entropy changes were derived using the temperature dependence of the free energy based on MD simulations at five different temperatures (see Methods Section). They show that complete pocket dehydration, starting from the average water occupancy, is entropically unfavorable. This observation agrees with the fact that solvent fluctuations are eliminated upon dehydration. At the same time, the free energy plateau reveals a perfect entropy−enthalpy compensation, indicating the favorable enthalpic effect due to the displacement of water molecules from the hydrophobic cavity interior to the bulk.

Cavity−Ligand Association

The obtained free energy profile for the ligand approaching the cavity is flat until ξ ∼ 0.4 nm. At shorter distances, it decreases monotonically indicating a steady mean force toward binding (Figure 4A). The free energy minimum is reached at ξ = −0.35 nm and corresponds to a direct contact between the ligand and the pocket bottom. The overall free-energy change upon binding, ΔG = −16.5 ± 0.6 kJ mol−1, is dominated by water-related contributions, ΔGW = −10.6 ± 0.6 kJ mol−1, as direct cavity−ligand interaction energy ΔUCL is only −5.9 kJ mol−1 (Figure 4B). Differently from what is found for other processes driven by hydrophobic interactions, nonpolar cavity−ligand binding is driven here by a large change in enthalpy, ΔH = −29.1 ± 17.3 kJ mol−1, with an opposing entropic contribution, −TΔS = 12.6 ± 17.3 kJ mol−1.

A relatively simple shape of the free energy profile covers substantial, compensating changes in enthalpy and entropy components as the ligand moves toward the cavity. A local maximum in H of 15.2 ± 6.1 kJ mol−1, mirrored by a minimum in −TS of −14.6 ± 6.1 kJ mol−1, is observed already around ξ = 0.55 nm. A comparison between water density maps for ξ = 0.55 and 1.1 nm (Figure 2) indicates that the reason for the observed changes is the partial destruction of the ligands second hydration shell, leaving the first hydration shell exposed to the weakly hydrated pocket region. Accordingly, a major contribution to the increasing enthalpy arises from the change in water−water interaction energy (U), with marginal effects due to cavity−water (UCW) and ligand−water (ULW) interactions.

The trends of enthalpic and entropic components suddenly invert as the ligand moves closer to the pocket, leading to a substantial decrease in H until ξ = 0.3 nm, with a subsequent plateau, and the formation of a −TS maximum (Figure 4A). These changes correspond to almost complete dehydration of the region between the cavity and the ligand (Figures 2 and 3A). The sizable enthalpy change observed is due to favorable contribution from water−water interactions that dominate unfavorable energetic effects of cavity and ligand dehydration. It is worth noting that no simultaneous change in free energy occurs, in agreement with the free energy profile for cavity dehydration described in the previous section. In addition, the change in the entropic component of PMF for the ligand moving from ξ = 0.55 to 0.3 nm (34 ± 11 kJ mol−1) is comparable with the entropic effect of changing the cavity occupancy from the average value to N = 1 (33 ± 28 kJ mol−1). This agreement might be partially fortuitous and should be treated semiquantitatively because of a number of factors. First, the cavity dehydration upon ligand binding and dewetting due to spontaneous fluctuations are different processes. Second, at ξ = 0.3 nm, the dehydration is not fully complete. Third, large uncertainties are associated with −TS values. Yet, the observed correspondence seems to confirm the overall interpretation of the described thermodynamic effects.

Further ligand translocation into the pocket corresponds to the gradual withdrawal of a hydrophobic object from an aqueous environment. It is accompanied by a steady decrease in free energy, initially driven by favorable changes in entropy. At the same time, no apparent change in enthalpy takes place, due to compensation between the increase in ULW and the decrease in UCW. For ξ < −2.5 nm, any extra ligand dehydration is limited. Increasingly favorable cavity−water interactions together with increasingly strong direct cavity−ligand attraction lead to a favorable enthalpy change. The moderate increase in the entropy component around the equilibrium binding distance (∼3 kJ mol−1) can be explained considering ligand immobilization due to direct contact with the wall.

Discussion

Our results demonstrate the determinant role of pocket (de)hydration upon hydrophobic cavity−ligand binding on the underlying thermodynamic signature. The favorable, driving enthalpy change results from the release of water molecules from the hydrophobic environment to bulk water. An opposing entropic component arises due to the elimination of solvent fluctuations inside the pocket. Such contributions dominate over the effects of extracting the small nonpolar ligand out of solution. Thus, cavity−ligand binding is enthalpy, rather than entropy, driven in our system and differs in this respect from other cases of nanoscopic hydrophobic association.[23,28,33] The hydrophobic association investigated herein is peculiar in that partial dewetting of one of the interacting partners occurs natively, already prior to binding due to the concave cavity geometry.

Remarkably, pocket dehydration displays virtually perfect enthalpy−entropy compensation, i.e., no net contribution to free energy. This is not expected to be a general phenomenon, as the free energy penalty upon dewetting depends on the specific enclosure size and the degree of hydrophobicity.[26,36,54−57] Larger or more hydrophilic pockets should be preferentially hydrated, while smaller hydrophobic pockets should remain natively dry.

Is our model system unique? A comprehensive study by Young and co-workers(34) focused on 14 protein structures selected from the Protein Data Bank (PDB) database for their large hydrophobic binding cavities. In explicit solvent molecular dynamics (MD) simulations drying was found in six cases, sometimes with few (usually two to five) water molecules fluctuating inside the cavity. Most of the reported cavities were long and narrow, as opposed to the hemispherical widely open pocket investigated herein. The driving forces for their hydration remain unknown, leaving open the interesting problem whether the architecture of the cavity may determine qualitative differences in thermodynamic signature for its interaction with water. A comparison of our results with the entropically hampered hydration found in spherical, closed cavities(36) indicates that the actual topography of the confining potential well and the degree of its openness toward bulk may be important in this respect. On the other hand, in line with our study, entropy-driven hydration was reported for a closed, nonpolar cavity of the I76A mutant of barnase.(41)

The (thermo)dynamics of cavity−ligand association addressed herein may depend on even a broader variety of subtle, balancing effects. Recent experimental studies shed light on nonpolar ligand binding to the poorly solvated pocket of the mouse major urinary protein-1 (MUP-1).[42−44] Despite the apparent hydrophobic character of the binding partners, association appears to be enthalpy-driven and accompanied by an unfavorable entropy change. Following a careful analysis, the authors attribute most of the enthalpic contribution to the effect of water displacement,(44) and the entropic penalty exclusively to the restriction of ligand and receptor degrees of freedom.(43) At the same time, in line with the argument by Dunitz, they assume that the release of bound water molecules is entropically favorable. Interestingly, however, the considered binding site seems to be hydrated by rather disorganized water.(44) This finding rises an intriguing question: is it legitimate to attribute the observed thermodynamic signature to loss in water entropy? As indicated by our study, such a scenario − in which cavity water is more entropic than bulk water − is possible.

The investigation of model cavity−ligand recognition to systems with varying physicochemical properties is being undertaken by means of the approach presented herein.(58)