Hydrophobic Association 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. Is this conventional picture directly transferable to biologi- cally relevant binding scenarios? They typically involve a receptor cavity, whose degree of hydration may vary signiﬁcantly 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 disordered39,40 with their displacement accompanied instead by an entropy loss,40 in line with a recent MD study.41 These important ﬁndings 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 ﬁve different temperatures. From the temperature dependence of the free energy, we obtain a complete picture of thermody- namic signature proﬁles 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 ﬁrst 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. J. Chem. Theory Comput., Vol. 6, No. 9, 2010 2867 2. Methods 2.1. Molecular Model and MD Simulations. The hemi- spherical cavity of 0.8 nm radius was embedded in a rectangular parafﬁn-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 ﬂat 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 nm46), and the TIP4P model was used for water molecules.47 Pairwise LJ interac- tions 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 ﬁlled with 1030 water molecules. Periodic bound- ary 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 bar51 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 2.2. Free Energy, Entropy, and Enthalpy Changes and Their Uncertainties. The Gibbs free energy for the system with the ligand at a given value reads G( ) ) W( ) + G(∞) (1) 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( ).