Cavity-ligand binding in a simple two-dimensional water model

By means of Monte Carlo computer simulations in the isothermal-isobaric ensemble, we investigated the interaction of a hydrophobic ligand with the hydrophobic surfaces of various curvatures (planar, convex and concave). A simple two-dimensional model of water, hydrophobic ligand and surface was used. Hydration/dehidration phenomena concerning water molecules confined close to the molecular surface were investigated. A notable dewetting of the hydrophobic surfaces was observed together with the reorientation of the water molecules close to the surface. The hydrogen bonding network was formed to accommodate cavities next to the surfaces as well as beyond the first hydration shell. The effects were most strongly pronounced in the case of concave surfaces having large curvature. This simplified model can be further used to evaluate the thermodynamic fingerprint of the docking of hydrophobic ligands.


Introduction
Due to the important influence of the confinement on the structural, thermodynamical, and phase equilibrium properties of simple fluids, these systems have been extensively studied recently [1][2][3]. The topic is also of great interest in the studies of biological systems in which many unanswered questions are concerned with hydration and dehydration of molecules of arbitrary shapes where water as a solvent finds itself confined close to the molecular surface.
One of the important issues that have begun to emerge is the effect of the shape, curvature, and roughness of a surface on its interaction with a ligand [4][5][6][7][8][9][10]. In order to understand the phenomena such as binding of the ligand to the receptor, or drugs to proteins, it is crucial to know the potential of the mean force between two biomolecules in question. Virtually all the binding sites in biology have a concave shape which imposes very particular geometrical constraints to the solvated water [11] and, therefore, the findings referring to the potential of the mean force between two spherical surfaces cannot be generalized to these problems.
To summarize the recent molecular dynamics studies on model systems of purely hydrophobic cavities [10,[12][13][14], it has been shown that water appears to be an active component in cavity-ligand association. An important impact of changes in water structure during the binding process has been noticed [10,14]. While the cavity hydration can be altered by changing its radius, in weakly hydrated cavity regions the reorganization of water molecules and suppression of the solvent fluctuations leads to enthalpy driven association [15]. A different interpretation of the computer simulation results have been given by Graziano [17]. His analysis shows that the Gibbs free energy gain upon association of a ligand in a concave hydrophobic cavity is mainly due to the decrease in the solvent-excluded volume, that translates in a gain of configurational-translational entropy of water molecules. This entropic driving force is masked by a large enthalpy gain associated with the reorganization of water-water hydrogen bonds upon association of the two nonpolar objects [17].
In this work we have focused on the investigation of the potential of the mean force between a hydrophobic ligand and a planar, concave, and convex purely hydrophobic surface with different curvatures. Our main interest was in the interpretation of the potential of the mean force through the water microstructural changes due to the confinement. To better visualize these effects we have chosen a simple two-dimensional Mercedes-Benz (MB) water model [16].
The paper is organized as follows: After this short introduction, the model and method are described. Next, the results are presented and discussed, and the conclusions are given in the end.

Model and method
The Mercedes-Benz model [16] was used to describe water molecules. In this model, the water molecule is represented as a two-dimensional Lennard-Jones (LJ) disk with three equally separated hydrogen bonding arms, and interacts with another water molecule through the potential U ww : (2.1) X i denotes the vector with coordinates and orientation of the i -th particle, and r i j is the distance between the centres of molecules i and j . Lennard-Jones 12-6 potential is: where LJ represents the depth of the potential well, and σ LJ is the Lennard-Jones diameter.
The hydrogen bond strength is a Gaussian function of the intermolecular distance and the angle between two hydrogen bonding arms: where G(x) is an unnormalized Gaussian function: The unit vectorî k represents the k-th arm of the i -th molecule (k = 1, 2, 3). Unit vectorû i j is the direction vector of the line joining the centres of the i -th and j -th molecule. Parameters HB = −1 and r HB = 1 determine the energy and the length of the optimal hydrogen bond. Values for the model parameters are as follows: σ LJ = 0.7 r HB = 0.7, LJ = 0.1 | HB | = 0.1, and σ = 0.085.
The hydrophobic ligand was represented as a Lennard-Jones disk of the same size (σ LJ = 0.7) as the water molecule.
To model a hydrophobic surface of a given curvature, a molecular wall was formed from Lennard-Jones discs (see figure 1). One surface was planar, while the others were bent. The interaction of the test water or hydrophobic ligand with the surfaces having positive curvatures (concave), and with negative curvatures (convex) was tested. The two radii used for concave surfaces were R * = 4.55 and 1.40, while for the convex surfaces the radii were 2.80 and 5.95. The hydrophobic particle and the water molecules interacted with the particles forming the hydrophobic surface through the Lennard-Jones potential [equation (2.2)]. The LJ parameters (σ LJ and LJ ) were the same for all particles in question.
Isobaric-isothermal Monte Carlo computer simulations [18] were performed for systems with 120 MB water molecules and a single hydrophobic ligand (test particle). The hydrophobic surface was placed in the centre of the simulation box. Reduced temperature, T * , and pressure, p * , of the systems were 0.20 and 0.19, respectively. Initial configuration of water molecules was chosen randomly and then equili- which started from equilibrated configuration. In one cycle, either rotation or displacement (probabilities for rotation and displacement being equal) of a water molecule was attempted. Maximal displacement and rotation were adjusted throughout the simulation, so that approximately one half of Monte Carlo moves were accepted. After every 120th cycle, an attempt to change the volume of the simulation box was made, where the maximal change was adjusted in the same way as maximal displacement/rotation. The pair correlation and angular distribution functions of water next to the surface were calculated using the histogram method, while the potential of the mean force between the hydrophobic ligand and the surface was calculated using the Widom's insertion technique [19].

Results and discussion
Here, we present the results of the Monte Carlo simulations of interaction of the hydrophobic ligand with the hydrophobic surface in water. All the results are given for T * = 0.2 and p * = 0.19. Figure 2 (a) shows the potentials of the mean force (PMFs) between the hydrophobic particle and the central particle  of the surface for convex surfaces, while the PMFs for concave surfaces are shown in figure 2 (b). The PMF with planar surface is shown for comparison in both panels. An important difference from the PMF between two hydrophobic discs in water (figure 1 of reference [6]) is observed. There is a larger difference in the values of the peaks belonging to the contact minimum (CM) and solvent separated minimum (SSM) compared to the homogeneous system (two hydrophobic discs); CM state is here much more stable than the SSM state. This suggests that the hydrophobic particle would prefer to stay close to the surface. Not much difference between the PMFs for different surface curvatures is observed in the convex cases, while in the case of strongly concave surface, the contact state is additionally stabilized and the position of the SSM is shifted further away from the surface.

13004-3
To interpret these results in view of water microstructure, we analysed the hydration of the surfaces more in detail. Characteristic simulation snapshots are shown in figure 3. A dehydration of the hydrophobic surface is noticed for all surfaces studied. The water molecules are pushed away and cavities are formed at the surface contact as well as beyond the first layer of water molecules. This is further confirmed by the water-surface radial distribution functions shown in figure 4. All g (r ) show a layered structure of water molecules away from the surface. There is no significant qualitative difference in the g (r ) between the planar surface case and convex bent surfaces. The value of the contact peak slightly decreases with an increasing curvature [see panel (a) of figure 4].
On the other hand, a qualitatively different behaviour is observed in the case of strongly concave surface. Here, water is completely pushed away from the cavity, forming the first hydration layer approximately two hydrogen bonds away from the surface [ figure 4 (b)]. This is in agreement with the pre-  vious studies by Baron et al. and Graziano [10,12,13,15,17]. Due to the incapability of forming hydrogen bonds, a water molecule in such a state is entropically stabilized. Further, the orientation of the surface water molecules belonging to the first hydration shell was examined. Figure 5 (a) shows the results for the convex surfaces, and panel (b) for the concave surfaces. The case of planar wall is shown in both panels for comparison. As already noticed by Southall et al. [20], a model water molecule close to the planar surface can no longer form all three hydrogen bonds with neighbouring water molecules. To compensate for this loss, it "wastes" one of the hydrogen bonds by pointing it towards the surface. This is seen from the angular distribution where the most probable orientation of the water molecule at the planar surface is the one pointing the hydrogen bond towards the surface (ϕ = 0°). The two hydrogen bonds pointing away from the surface participate in the formation of cavities facilitating a favourable SSM in the PMFs.
As the surface adopts the convex shape, the angular preference of the water to the surface is much less expressed. For the case in panel (a) of figure 5, we can conclude that angular preferences of water have a negligible contribution to the shape of the PMFs.
The situation gets reversed in the case of strongly concave surface [triangles in figure 5 (b)]. The preferential angle of water molecule in this case is the one where the hydrogen bonding arm points at an angle 60°with respect to the centre-centre coordinate [see the insert in panel (a)], suggesting the strongly obstructed hydrogen bond formation due to the lack of space in the cavity. This way a water molecule was stabilized by an increase of the rotational entropy.

Conclusions
The results for the potential of the mean force between hydrophobic solute and hydrophobic surface of various shapes in the model water solutions were presented. The results were analysed in view of hydration water structure and orientation. All the hydrophobic surfaces studied showed dewetting which was most pronounced for the bent concave surface. The angular orientation of the water molecules facilitated the formation of cavities which stabilized the non-covalent direct and water separated binding of the hydrophobic molecule to the surface. In the future work, temperature dependence of the phenomena will be investigated and the results will be used for the enthalpy-entropy decomposition of the free energy of ligand-surface interaction.