Dewetting-controlled binding of ligands to hydrophobic pockets

We report on a combined atomistic molecular dynamics simulation and implicit solvent analysis of a generic hydrophobic pocket-ligand (host-guest) system. The approaching ligand induces complex wetting/dewetting transitions in the weakly solvated pocket. The transitions lead to bimodal solvent fluctuations which govern magnitude and range of the pocket-ligand attraction. A recently developed implicit water model, based on the minimization of a geometric functional, captures the sensitive aqueous interface response to the concave-convex pocket-ligand configuration semi-quantitatively.

The water-mediated interaction between a ligand and a hydrophobic binding pocket plays a key role in biomolecular assembly processes, such as protein-ligand recognition [1,2,3,4,5,6], the binding of the human immunodeficiency virus (HIV) [7] or the dengue virus [8] to human cells, the inhibition of influenza virus infectivity [9], or in synthetic host-guest systems [10]. Experiments and explicit-water molecular dynamics (MD) simulations suggest that the concave nature of the host geometry imposes a strong hydrophobic constraint and can lead to very weakly hydrated pockets [1,2,3,4,5,6,11,12], prone to nanoscale capillary evaporation triggered by the approaching ligand [4,5,13]. This so called 'dewetting' transition has been also observed in other (protein) geometries [11,12].
It has been speculated that dewetting may lead to a fast host-guest recognition accelerating the hydrophobic collapse and binding rates of the ligand into its pocket [1,4,5]. A deeper physical understanding of these sensitive hydration effects in hydrophobic recognition is still elusive.
On the coarse-grained modeling side, the thermodynamics of molecular recognition is typically approached by surface area (SA) models [14]. A major flaw of these implicit solvent models is that the aqueous interface around the macromolecules is predefined (typically by rolling a probe sphere over the van der Waals surface) and is therefore a rigid object that cannot adjust to local energetic potentials and changes in spatial molecular arrangements.
In particular, the dewetting transition, which is highly sensitive to local dispersion, electrostatics, and geometry [11,12,15], can, per definitionem, not be captured by SA type of models. Their qualitative deficiency to describe the hydrophobic pocket-ligand interaction in proteins [16], pocket models [13], or dewetting in protein folding [17] is therefore not surprising.
In this letter, we combine explicit-water MD simulations and the variational implicit solvent model (VISM) [15] applied to a generic pocket-ligand model. The simulations show that the approaching ligand first slightly stabilizes the wet state in the weakly hydrated pocket, whereas, upon further approach, bimodal fluctuations in the water occupancy of the pocket are induced, followed by a complete pocket dewetting. The onset of fluctuations defines the critical range of pocket-ligand attraction. The VISM calculation, based on the minimization of a geometric functional, reproduces the bimodal hydration and explains it by the existence of distinct metastable states which correspond to topologically different water interfaces. As opposed to previous (SA type of) implicit models, VISM captures the range of the pocket-ligand interaction semi-quantitatively. Strikingly, the observed nanoscale phenomena can be thus explained by geometric capillary effects, well-known on macroscales [18]. Explicit inclusion of dispersion interactions and curvature corrections, however, seem to be essential for an accurate description on nanoscales.
Our generic pocket-ligand model consists of a hemispherical pocket embedded in a rectangular wall composed of neutral Lennard-Jones (LJ) spheres interacting with U LJ (r) = 4ǫ[(σ/r) 12 − (σ/r) 6 ]. The atoms are aligned in a hexagonal closed packed arrangement with a lattice constant of 1.25Å. The LJ parameters are chosen to model a paraffin-like material and are ǫ p = 0.03933 kJ/mol and σ p = 4.1Å [13,27]. We consider two different pocket radii: R = 5 and 8Å, which we refer to in the following as 'R5' and 'R8' systems. The ligand is taken as a methane (Me) represented by a neutral LJ sphere with parameters ǫ = 0.294 kJ/mol and σ = 3.730Å. It is placed at a fixed distance d from the flat part of the wall surface (z = 0), along the pocket symmetry axis in z-direction, see Fig. 1 (a) for an illustration. The explicit-water MD simulations are carried out with the CHARMM package [19] employing the TIP4P water model in the NV T ensemble, two anti-symmetric walls with thickness of 7.5Å in a surface-to-surface distance of 30Å in a rectangular box with lengths L x = L y ≃ 34Å and L z = 100Å, and 3D particle mesh Ewald summation.
In equilibrating simulations, the volume V of the system was varied until the density in the center of the slab matched the bulk density of TIP4P water at a pressure of P = 1 bar and temperature T = 298 K. More technical details of the simulation setup can be found elsewhere [13,27]. A MD simulation snapshot is shown in Fig. 1 The VISM was introduced in detail previously [15] and applied to the solvation of nonpolar solutes [20]. Briefly, let us define a subregion V void of solvent in total space W, for which we assign a volume exclusion function v( r) = 0 for r ∈ V and v( r) = 1 else. The volume V and interface area S of V can then be expressed as functionals and S[v] = W d 3 r |∇v( r)| = ∂W dS, and the solvent density is ρ( r) = ρ 0 v( r), where ρ 0 is the bulk value. The solvation free energy G is defined as a functional of the geometry v( r) of the form [15] where γ lv is the liquid-vapor interface tension, δ the coefficient for the curvature correction of γ lv in mean curvature H( r), and U( r) = Ns i U i LJ ( r − r i ) sums over the LJ interactions of all N s solute atoms at r i (ligand+wall atoms) with the water. The δ-term in (1) has been used in scaled-particle-theory [21] for convex solutes only, generalized capillary theory [23], and in the morphometric approach applied to the solvation of model proteins [22]. The minimization δG[v]/δv = 0 leads to the partial differential equation (PDE) [15] which is a generalized Laplace equation of classical capillarity [18,23] extrapolated to microscales by dispersion and the local Gaussian curvature K( r). The PDE (2) is solved using the level-set method which relaxes the functional (1) by evolving a 2D interface in 3D space and robustly describes topological changes, such as volume fusion or break-ups [20,24]. The free parameters chosen to match the MD simulation are P = 1 bar, γ lv = 59 mJ/m 2 for TIP4P water [25], and ρ 0 = 0.033Å −3 . The coefficient δ is typically estimated to be between 0.8 and 1Å for various water models around convex geometries [11,26], while VISM was able to predict well the solvation free energies of simple solutes for δ = 1Å [20] which we use in the following.
We consider ligand positions from d = 11Å to the distance of nearest approach to the pocket bottom. The latter is defined as corresponding to a wall-ligand interaction energy of 1 k B T and is d ≃ −1.8 and -3.8Å for the R5 and R8 system, respectively. We define the water occupancy N w of the pocket by the number of oxygens whose LJ centers are located at z < 0. Considering the probability distribution P (N w ), we obtain the free energy as a function of N w by G(N w ) = −k B T ln P (N w ) + G ′ . Without the ligand (effectively for d 9Å), the MD simulation reveals that the R5 pocket is in a stable dry state with occupancy N w ≃ N dry = 0, despite the fact that a few water molecules fit in and consistent with experiments on an equally sized protein pocket [6]. The R8 system, however, is found to be weakly hydrated. The G(N w ) distribution shown for d = 9Å in Fig. 2 reveals an almost barrierless transition between wet and dry states. The metastable wet state comprises N w ≃ 9 = N wet water molecules which roughly corresponds to bulk density.
The approaching ligand considerably changes the G(N w ) distribution in the R8 system.
As plotted in Fig. 2, for d = 6.5Å the free energy exhibits a minimum at the wet state which is slightly stabilized (by ≃ 0.4 k B T ) over the dry state. The function G(N w ) develops, however, concave curvature for N w ≃ 0 indicative of the onset of a thermodynamic instability. Indeed, upon further approach of the ligand (d = 5.5Å) a local minimum forms at the dry state. It becomes a stable, global minimum at the critical distance d c ≃ 4.5Å.
The now metastable wet state completely vanishes for d 0Å, where we find a free energy difference between the wet and dry state of G(N dry ) − G(N wet ) ≃ 5k B T . By investigating the water density distribution (Fig. 2, right panel)  and 2s-dry, indicating a very stable dry pocket in agreement with the MD simulations and experiments [6]. These results demonstrate that VISM captures the dewetting transition, and the final interface is relaxed into (meta)stable states representing (local) free energy minima. This is in physical agreement with the bimodal behavior observed in the MD simulation and is further quantified in the following.
The minimum VISM free energy (1) vs. d is is plotted for R8 in Fig. 4 Changing the curvature parameter δ shows that this failure can be attributed to a too high energy penalty for concave interface curvature (a too large δ for H < 0) which favors pocket dewetting. It thus appears that the simple curvature correction applied breaks down and is In summary, the geometry-based VISM is the first implicit solvent model that captures the multi-state hydration observed in simulations and experiments [12] and highlights the significance of interfacial fluctuations [29] in hydrophobic confinement where the free energy can be polymodal. Pocket dewetting may be regarded as the rate-limiting step for proteinligand binding as found in folding [11].

The wall-water interaction
In order to construct hydrophobic walls we considered a paraffine-like material of 0.8 g/cm 3 density composed of CH 2 units. Assuming a hexagonally closed-packed arrangement, the given density requires a lattice constant of 3.5Å which is too coarse to produce a relatively smooth hemispherical pocket. Thus, we reduced the lattice constant to 1.25Å while at the same time adjusting the Lennard-Jones potential parameters of the wall pseudo-atoms to reproduce the original paraffine wallwater interaction energy (see inset to Fig. 1) that was obtained with the united atom OPLS parameters for CH 2 units [1]. The wall-water interaction energy was calculated by simply averaging the interaction energy over planes at constant z. The height of the first peak in the wall-water density profile from our MD simulations (Fig. 1) is within the range (1.3 to 1.6) observed in all atom MD simulations of hydrocarbon-water interfaces [1,2,3], suggesting that the walls indeed closely resemble a paraffine-like material.
where N [v] = 0 for dry-type solutions and N [v] = N wet = 9 for wet-type solutions. A comparison with MD results (Fig. 2) shows the correct qualitative behavior with a transient increase in pocket wetting at d ≃ 6Å, followed by ligand induced dewetting. The quantitative discrepancy is probably due to a) the approximations in the VISM functional leading to over-stabilization of dry state relative to wet state, and b) including only local G[v] minima in the ensemble average while omitting intermediate states of not much higher free energy. We expect the qualitative agreement to improve upon inclusion of interface thermal fluctuations into VISM.