1. Introduction
Protein maintains its structure and function upon solvation in water. There is increasing evidence supporting that water not only acts as a solvent but also actively participates in many biological processes [1, 2]. For example, it has been shown that the protein hydration water has a significant impact on protein dynamics [3–5], protein-ligand binding [6, 7], protein stability [8–11] and the catalytic efficiency of enzymes [12]. Therefore, the solvation of protein is key to understanding the biological function of proteins.
Water being able to form a connected hydrogen-bond (H-bond) network with locally favored tetrahedral symmetry is the most unique and anomalous solvent in nature [13–16]. It has been shown that the local tetrahedral ordering is responsible for both the thermodynamic and dynamic anomalous behaviors of water [17–20]. The presence of solute inevitably perturbs the tetrahedral structure of water. In 1959 Kauzmann proposed that the water structural ordering around hydrophobic solutes is the origin of the hydrophobic interaction which serves as the key driving force of the protein folding and aggregating in aqueous solutions [21]. It's now well accepted that the hydrophobic interaction is entropic in its origin, but how to characterize the underlying water structure in the vicinity of proteins has remained a major challenge so far.
The structure of protein hydration water has been intensively studied by using various structural descriptors. However, either experiments or simulations report contradictory effects of protein on the hydration water structure. For example, Shen et al reported that hydration water of most amino acids has higher tetrahedral order than bulk water by using Raman multivariate curve resolution spectroscopy [22]. Enhanced H-bonding structure of protein hydration water has also been found by Fourier transform infrared spectroscopy [23] and femtosecond surface sum frequency generation spectroscopy [24]. These results are supported by molecular dynamics (MD) simulations that reported a significantly structured hydration water layer around a lysozyme protein [25, 26].
In contrast, neutron Brillouin measurements combined with MD simulations reported that lysozyme protein breaks the tetrahedral order of hydration water [27]. X-ray scattering experiments [28] and MD simulations [29, 30] also detected reduced tetrahedral order in the protein hydration layer. Moreover, Merzel and Smith found that the hydration water of lysozyme is 15% denser than bulk water [31]. Since density is anticorrelated to the local structural order of water [32], the increased density supports the depletion of tetrahedral order in the protein hydration layer.
The effect of protein on the local structural ordering of hydration water has remained elusive so far [2]. The difficulty arises from the fact that neither the translational nor the rotational symmetry preserves at the protein/water interface, and thus, traditional structural descriptors targeting the tetrahedral order may not be suited for protein hydration water at the interface [25]. In this work, we analyzed the water H-bond structure and applied a newly developed structural descriptor to protein hydration water. We find that the structural characterization focusing on the H-bond network unambiguously detects enhanced local structural ordering of the protein hydration water. This work not only opens a new door to the structural characterization of protein hydration water but also provides microscopic evidence supporting Kauzmann's seminal idea on the hydrophobic interaction.
2. Methods
In this study, we take the hen egg white lysozyme as the model protein, since it has been widely studied as an archetype protein in both experiments and simulations. The lysozyme protein contains 129 residues and the initial structure is obtained from the protein data bank (ID: 1IEE) [33]. The CHARMM36 force field [34] was adopted to describe the interactions of protein and the water was modeled by the TIP4P/2005 model [35]. A lysozyme protein was solvated in a cubic box of 63 082 water molecules. Eight chloride ions were added to keep the charge neutrality of the system. The box is around 124 × 124 × 124 Å with the periodic boundary condition applied in all directions. The system was equilibrated at 300 K and 1 bar for 5 ns and followed by another 2.1 ns NVT equilibration run at 300 K with the volume determined from the NPT run. Then a production run was performed in an NVT ensemble at 300 K for 2.4 ns and the configurations were sampled every 0.2 ps. All the bonds with hydrogen atoms were constrained by the LINCS algorithm. A time step of 2 fs was adopted for the simulations. The temperature and pressure were kept constant by using the Nose–Hoover thermostat and the Parrinello–Rahman barostat, respectively. The van der Waals and the electrostatic interactions in the real space were truncated at 12 Å and the electrostatic interactions in the reciprocal space were treated by the fast smooth particle-mesh Ewald method. Simulation of pure water was carried out in a system of 27 000 TIP4P/2005 water molecules at 300 K for 2.4 ns with the other parameters the same as the protein simulations. All the simulations were performed by using the GROMACS (2019.4) package [36, 37].
3. Results and discussion
In pure water, molecules favor the tetrahedral arrangement of neighboring molecules. The degree of the tetrahedral order can be described by a parameter q as [38, 39]
$\begin{eqnarray}q=1-\displaystyle \frac{3}{8}\sum _{i=1}^{3}\sum _{j=i+1}^{4}{\left(\cos {\theta }_{{ij}}+\displaystyle \frac{1}{3}\right)}^{2},\end{eqnarray}$
where θij is the angle formed by two vectors connecting the central molecule and its nearest neighbors i and j, and the summation runs over all the combinations of the four nearest neighbors. It takes a value of 0 and 1 for a random and a perfect tetrahedral configuration, respectively. The tetrahedral parameter q has been widely used to characterize the tetrahedral order of protein hydration water [27, 29, 30]. In pure water, the parameter q is defined by using the oxygen atoms only. However, as pointed out by Accordino et al, the parameter q may not be suited for water at the interface, because interfacial water may not often have four neighbors in the first coordination shell [25]. The protein N and O atoms that are able to form H-bonds with water are often involved in the definition of q to compensate for the loss of neighboring water molecules at the protein/water interface. The tetrahedral parameter q targeting the rotational (tetrahedral) symmetry is determined solely by the angular distribution of neighbors. However, the presence of protein inevitably breaks the translational and rotational symmetry of the water arrangement at the interface.There are many other structural descriptors focusing on the translational order of water, such as d5 [40] and the local-structure index [41]. These structural descriptors have been successfully applied to the characterization of the translational order of pure water. However, none of the above-mentioned structural descriptors, including the tetrahedral parameter q, consider the H-bond formation in their definitions. Since H-bond formation is the essential driving force for water structuring, characterization of the H-bond network in the protein hydration layer is crucial to reconcile the discrepancy in the structural description of protein hydration water.
To find the protein hydration water, we adopted the cutoff method by which a water molecule is selected as protein hydration water if it is within 5 Å of at least one carbon atom of the protein [42]. Persson and coworkers have shown that this cutoff method is able to detect protein hydration water efficiently and accurately [42]. Figure 1(a) shows the distribution of coordination number nFS (the number of water molecules in the first coordination shell) of protein hydration water and bulk water. We can see that the protein hydration water has a similar distribution shape compared to bulk water, but the peak position shifts from nFS = 5 for bulk water to nFS = 4 for protein hydration water. This shift of the distribution corresponds to the reduction of the water coordination number from 5.09 for bulk water to 4.18 for protein hydration water (table 1), which is ascribed to the confinement effect induced by the presence of protein. Figure 1(b) displays the distribution of the number nHB of H-bonded water neighbors per water molecule for protein hydration water and bulk water. Here, two water molecules are considered as H-bonded if their oxygen–oxygen distance is smaller than 3.5 Å, and the H–O⋯O angle is less than 30° [43, 44]. The presence of protein not only shifts the peak position of the distribution from nHB = 4 for bulk water to nHB = 3 for protein hydration water but also changes the shape of the distribution. This can also be seen from table 1 that each protein hydration water loses 0.58 H-bonds, compared to a loss of 0.91 neighbors on average, in the presence of protein, which suggests that the effect of protein is not only spatial confinement but also leads to the reorganization of water H-bond network.
Figure 1. The distribution of (a) the coordination number nFS and (b) the number nHB of H-bonded water neighbors per water molecule for protein hydration water (black squares) and bulk water (blue circles). |
Table 1. The average value of the structural descriptor ζ, the coordination number nFS, the number of H-bonded water neighbors nHB, the number of non-H-bonded water neighbors Δn, the H-bond length rHB, and the H-bond strength εHB of protein hydration water and bulk water obtained from our simulations. The standard deviations of the structural descriptors are shown in the parentheses. |
| Water type | ζ (Å) | nFS | nHB | Δn | rHB (Å) | εHB (kJ mol−1) |
|---|---|---|---|---|---|---|
| Protein hydration water | 0.37 (0.40) | 4.18 (1.25) | 3.08 (0.94) | 1.10 (1.12) | 2.86 (0.17) | 18.91 (5.68) |
| Bulk water | 0.26 (0.35) | 5.09 (1.07) | 3.66 (0.73) | 1.43 (1.25) | 2.87 (0.18) | 18.92 (5.66) |
To characterize the effect of protein on water H-bond structure, we calculated the number of non-H-bonded water neighbors in water's first coordination shell which is defined by the following relation,
$\begin{eqnarray}{\rm{\Delta }}n={n}_{\mathrm{FS}}-{n}_{\mathrm{HB}}.\end{eqnarray}$
The distributions of Δn for protein hydration water and bulk water are shown in figure 2. As we can see, protein significantly promotes the formation of a fully H-bonded first coordination shell (Δn = 0) by ∼40% and depletes the coordination shell with non-H-bonded neighbors, compared to bulk water. This result indicates that protein promotes the water H-bond structure in its hydration shell. Sciortino et al have demonstrated that the presence of non-H-bonded molecules, which may be treated as ‘defects' in the first coordination shell, effectively enhances molecular mobility in liquid water [45]. Thus, the promotion of the fully H-bonded first coordination shell should slow down the mobility of protein hydration water, which agrees with previous simulation and experimental results [46, 47].
Figure 2. The distribution of the number Δn of non-H-bonded water neighbors in water's first coordination shell per water molecule for protein hydration water (black squares) and bulk water (blue circles). |
Recently, Russo and Tanaka proposed a new structural descriptor ζ to characterize the local translational order of liquid water [48]. The descriptor ζ measuring the depth of non-H-bonded water penetrating into the first coordination shell is defined for each water molecule as
$\begin{eqnarray}\zeta ={d}_{\mathrm{nhb}}-{d}_{\mathrm{hb}},\end{eqnarray}$
where dnhb and dhb are the distance from the closest non-H-bonded water and the distance from the furthest H-bonded water to the central molecule, respectively. A small ζ around 0 corresponds to a disordered structure with a penetrated non-H-bonded molecule in the first coordination shell, whereas a relatively large ζ suggests a translationally ordered water structure with a fully H-bonded first coordination shell. Incorporating the H-bond information, the ζ parameter has been successfully exploited to characterize the local structural ordering in pure water [48–50]. Figure 3 plots the distribution of ζ for protein hydration water and bulk water. Clearly, the protein hydration water has a broader distribution than bulk water, which may be attributed to the topological and chemical heterogeneities of the protein surface. Moreover, the ζ distribution for protein hydration water shifts towards a large ζ value, compared to bulk water. Accordingly, the average ζ value increases by 42% from 0.26 Å for bulk water to 0.37 Å for protein hydration water (table 1). The ζ distribution clearly demonstrates that the protein hydration water is structurally more ordered than bulk water, in agreement with the above analysis of the H-bond network (figures 1 and 2).
Figure 3. The distribution of the structural descriptor ζ for protein hydration water (black solid line) and bulk water (blue dash line). |
Besides H-bond structure, we also investigated the effect of protein on the length rHB and strength εHB of individual H-bond. Here, rHB and εHB are defined as the oxygen–oxygen distance and the interaction energy (in absolute value) of two H-bonded water molecules. Figures 4(a) and (b) show the distribution of rHB and εHB, respectively. In contrast to the significant impact of protein on water H-bond structure, the presence of protein turns out to have negligible influence on either the length or the strength of water–water H-bond statistically. This result suggests that the protein promotes water's local structural ordering through the reorganization of the H-bond network, rather than perturbing the strength of individual H-bond.
Figure 4. The distribution of (a) the H-bond length rHB and (b) the H-bond strength εHB for protein hydration water (black solid line) and bulk water (blue dash line). In panel (b) the two distributions overlap. |
The structure and dynamics of water in the near vicinity of the protein are rather heterogeneous [51, 52]. It has been shown that geometric topology [53, 54], charge distribution [55], chemical nature [56] and concentration [57] of the protein all affect the structure of protein hydration water. In figure 5 we show the snapshot of the lysozyme protein and its hydration layer. It can be seen that the cutoff method [42] accurately selected the hydration water at the protein/water interface. To illustrate the spatial heterogeneity of water's local structural order, we calculate the ζ parameter of each hydration water and show the instantaneous value of ζ by the color of water molecules in figure 5(b). As clearly indicated by the ζ parameter, the local structural ordering of protein hydration water takes place heterogeneously on the protein surface. Moreover, the protein hydration water molecules with similar ζ values tend to aggregate into small patches on the protein surface, suggesting that the local structural ordering is not random but takes place in a cooperative manner. We note that the value of ζ fluctuates with time due to the thermal fluctuations of water structure at finite temperatures. Therefore, we calculated the average value of ζ for each residue, ⟨ζ⟩, that is defined by averaging the ζ value over all the water molecules in contact with the residue and over time. Here, a water molecule is considered in contact with a residue if it is in the hydration shell of the protein and the residue is the closest one to that water molecule. The average value ⟨ζ⟩ provides a measure of the degree of water structuring in the vicinity of each residue. We plot the spatial distribution of ⟨ζ⟩ in figure 5(c). It can be seen that the structure of hydration water is indeed heterogeneous on the protein surface and this structural heterogeneity is strongly correlated with the residues on the protein surface. Understanding the origin of the spatial heterogeneity of water structuring on the protein surface and its link to the structure and chemical nature of the residues is of great interest for future study.
Figure 5. (a) The snapshot of lysozyme protein in water at 300 K. (b) The snapshot of hydration water (small balls) on the surface of lysozyme protein (grey cloud) at 300 K. The water molecules are colored by the instantaneous value of the structural descriptor ζ. (c) The surface of lysozyme protein at 300 K. The residues on the surface are colored by the value of ⟨ζ⟩ that is defined by averaging the ζ value over all the water molecules in contact with the residue and over time. The blue and red color represents ordered and disordered water, respectively. The color bars for ζ and ⟨ζ⟩ are shown in panels (b) and (c), respectively. The H-bonds formed between hydration water molecules are shown by black sticks. |
4. Summary
We have studied the effect of lysozyme protein on the structure of hydration water by all-atom MD simulations. Previous studies characterize the protein hydration water by structural descriptors focusing on either tetrahedral or translational order neglecting H-bond information. In this work, we have focused on the structure of water's H-bond network at the protein/water interface. We find that the protein facilitates the formation of a fully H-bonded first coordination shell of water in absence of any penetrating non-H-bonded molecules (‘defects') on the protein surface. Moreover, the presence of protein tends to deplete the disordered water structure with non-H-bonded molecules in water's first coordination shell. Applying a newly developed translational structural descriptor ζ that explicitly takes H-bond formation into account, we find that the presence of protein promotes the ζ value by 42% for the protein hydration water compared to bulk water. This result, together with the analysis of the H-bond network, clearly demonstrates the significant development of the local structural order of water at the protein/water interface. This work highlights the essential role of H-bonding in the structural characterization of the interfacial water [50] and provides clear microscopic evidence for the water structural ordering around the protein that underlies the essential hydrophobic interactions in biological systems [21]. As an archetype protein, lysozyme contains various kinds of residues (polar, non-polar, positively charged, and negatively charged ones) and forms different types of protein structures (α-helices, β-sheets, and loops). Thus, the results obtained from lysozyme protein are expected to be relevant in general for other proteins as well. This work provides new insights into the microscopic structural characterization of protein hydration water and is fundamental to the understanding of the solvation of biomolecules in water.