Revealing the Frank–Evans “Iceberg” Structures within the Solvation Layer around Hydrophobic Solutes

Using computer simulations, the structural properties of solvation water of three model hydrophobic molecules, methane and two fullerenes (C60 and C80), were studied. Systems were simulated at temperatures in the range of 250–298 K. By analyzing both the local ordering of the molecules of water in the solvation layers and the structure of hydrogen bond network, it is shown that in the solvation layer of hydrophobic molecules, ordered aggregates consisting of water molecules are formed. Even though it is difficult to define the exact structure of these aggregates, their existence alone is clearly noticeable. Moreover, these aggregates become more pronounced with the decrease of temperature. The existence of the ordered aggregates around the hydrophobic solutes complies with the concept of “icebergs” proposed by Frank and Evans.


Details of the preparation of the systems and the computer simulations
There were some common settings of all production runs for all simulated systems. The particle-mesh Ewald method was used for electrostatic interactions. The lengths of chemical bonds involving hydrogen atoms were fixed using SHAKE. A 1.2 nm cutoff for nonbonding interactions was used. Simulations were performed under NPT conditions. Constant temperature (298 K) was kept by the weak coupling to an external bath (τ T = 1.0 ps); constant pressure (1 bar) was kept by the weak coupling method (τ p = 1.0 ps). Details of the simulation protocol varied between the systems.
Simulation procedure for the system with the lysozyme molecule In order to ensure that the presence of ions is not affecting the results of the analysis of the solvation shell of the protein molecule, two separate simulation systems were prepared. In the first one, ions were distributed randomly in the system and no restrains were applied on their positions. In the second one, ions were placed close to the corners of the simulation box, to maximize their distance from the surface of the lysozyme molecule. The positions of the chloride anions were restrained at the reference points using harmonic potential, with the value of the force constant equal to 2.0 kJ/(mol · A 2 ). In order to prevent the protein molecule from moving towards the ions, the backbone carbon atoms were also restrained, with the use of force constant equal to 0.02 kJ/(mol · A 2 ). Both systems were simulated using the same conditions as in the other systems (see main text). The analysis showed that the placement of the chloride ions had no significant influence on the results ( Figure S1). Figure S1: The ordering map illustrating the temperature changes of the degree of structural order of the solvation water of the studied solutes: the basal plane of ice, fullerenes C60 and C80, methane, lysozyme and chloride ions. The color gradient indicates the direction of the temperature change (the lighter the color the lower the temperature). The results for restrained lysozyme with 0.5 nm-thick solvation layer (LZ), restrained lysozyme with 0.45 nm-thick solvation layer (LZ_0.45) and unrestrained lysozyme with 0.45 nm-solvation layer (LZ.ur_0.45) are shown for comparison.

Radial distribution function
The first step in the calculations was to determine the extent of the first solvation shells of the molecules chosen for the study. It was achieved on the basis of the "surface" density distribution function. The distance at which this function reaches a minimum is considered to be the boundary of the first solvation layer of the molecule of a solute. The function was determined as the ratio of the number of solvation water molecules, N solv , at a given distance r from the surface of the molecule of the solute, to the number of N bulk bulk water molecules at the same distance in a "reference system" in which the molecule under consideration was inserted into a box filled with bulk water: The course of the functions is presented in the figure below. The position of the first minima of the functions determines the extent of the first solvation layers -they are located at a distance of about 0.50 nm from the surface of the molecules of the hydrophobic solutes. The course of the function in the case of the lysozyme is less straightforward: the first maximum of the function g(r) clearly consist of two overlapping peaks. It reflects the complex nature of the interaction of water molecules with the protein surface. The concept of the order parameters -two-particle correlation function The concept of the order parameters starts with the two-particle correlation function. Two-particle correlation function is -by definition -a ratio of two probabilities: the probability of finding a second molecule in a certain position (distance, configuration and orientation) relative to the central molecule in a studied system (dw ) and the probability of finding a second molecule in the same position in a system where the positions of molecules are not correlated (dw).
The position of one water molecule relative to the other can be characterized by their distance (r) and five angles (φ, θ, α, β, γ), describing the distance and orientation of two coordinate systems. The coordinate systems is ascribed to a water molecule as depicted in S3. In our case, only the distance and configuration (described by r and the angles φ and θ) were considered (see also Figure S3).
The first probability can be calculated by dividing the number of molecules located in a certain space volume around the central molecule (dN (r, φ, θ)) by the total number of molecules in the studied system (N − 1 ≈ N ). The second probability is calculated by dividing the volume of the considered fragment of the space (dV (r, φ, θ)) by the total volume of the studied system (V ). Figure S3: Schematic representation of the variables needed for the calculation of a two-particle correlation function, when only distance and configuration are analyzed. Axes Ox, Oy, Oz are related to the central water molecule. The coordination system associated with each molecule is defined as shown in the right side of the figure.
Therefore, the two-particle correlation function can be calculated as follows: (1) Two-particle correlation function can then be decomposed into translational and configurational contribution: where dN (φ, θ|r) is a number of molecules located in the area described by the φ and θ angles (in relation to the central molecule), assuming that they are within the distance range of r, r + dr from the central molecule. The functions g (2) tra (r) and g (2) con (φ, θ|r) are then used to calculate the order parameters p tra and p con , as discussed in the main text.

Interpretation of the values of the order parameters
The values of the p tra parameter (and the p con parameter) depend on the integral from the two-particle correlation function g (2) (r, φ, θ). The values of the two-particle correlation function for water in solvation shell depend on how the space around some water molecule is occupied by other water molecules. Therefore, the values of this function for water in solvation shell can be different from the values for bulk water, if preferable distances and mutual orientations of water molecules differ. However, some space can be unoccupied by other water molecules because it is already occupied by the solute -this also influences the values of the p tra and p con parameters. We hope to diminish the influence of the second factor on the results by subtracting the results for bulk water filling a fictitious solvation shell (Kuffel, A.; Zielkiewicz, J. The Importance of the Shape of the Protein-Water Interface of a Kinesin Motor Domain for Dynamics of the Surface Atoms of the Protein. Phys. Chem. Chem. Phys. 2012, 14, 5561). We have to be aware that the geometry of the solvation layer (including the cutoff used for the construction of the artificial solvation layer filled with bulk water, here chosen as equal to 0.17 nm -as described in the main text) influences the results. It is the reason why we are more interested in comparing temperature changes of the order parameters (their extent and direction observed for various systems) rather than in analysing their absolute values.
The values of these parameters cannot be directly translated to any specific spatial structure. Solvation water molecules in solvation layer of the studied solutes are restrained in their own way, depending on the nature of the solvated surface, but it is possible that different structural arrangements of molecules can be characterized by similar values of the order parameters. For example, results presented on the ordering map show that the values of translational and configurational order parameters obtained for solvation water of fullerenes and ice are not that different at higher temperatures. At the same time there are evidence that the structure of solvation water is not the same in these cases (see Figures 2 and 4 in the main text). In our opinion, there is no contradiction in these results. Fullerenes and ice differ significantly in terms of the chemical composition and spatial structure of their surfaces. Fullerenes, because of their hydrophobic character, are unable to form hydrogen bonds with the molecules of the solvent -water S6 molecules in the solvation shell are therefore restrained to a limited range of mutual arrangements if they want to maximize their hydrogen bonding (what is achieved at the cost of some deformation of the tetrahedral arrangement of the hydrogen-bonded molecules - Figure 4). Near the ice block, on the other hand, water molecules are induced to adopt certain positions relatively to the ordered water molecules in the ice block. In both cases increased translational and configurational ordering (compared to bulk water) can be observed, even though its origin appears to be different.
The analysis of the water-water hydrogen bonds -the δ angles Figure S4: Schematic representation of the definition of the δ angle. Orange points indicate the location of oxygen atoms of the water molecules in the ring. Dark blue vector connects the centre of the ring and the geometrical centre of a set of ten atoms located on the surface of solute (or of all the atoms in the case of methane molecule). Light blue vector is a vector starting at the centre of the ring and normal to the ring plane.