Characterizing Hydration Properties Based on the Orientational Structure of Interfacial Water Molecules
CCharacterizing Hydration Properties Based on the OrientationalStructure of Interfacial Water Molecules
Sucheol Shin and Adam P. Willard ∗ Department of Chemistry, Massachusetts Institute of Technology,Cambridge, Massachusetts 02139, USA
Abstract
In this manuscript, we present a general computational method for characterizing the molecularstructure of liquid water interfaces as sampled from atomistic simulations. With this method, theinterfacial structure is quantified based on the statistical analysis of the orientational configurationsof interfacial water molecules. The method can be applied to generate position dependent maps ofthe hydration properties of heterogeneous surfaces. We present an application to the characteri-zation of surface hydrophobicity, which we use to analyze simulations of a hydrated protein. Wedemonstrate that this approach is capable of revealing microscopic details of the collective dynamicsof a protein hydration shell. a r X i v : . [ c ond - m a t . s o f t ] A ug n the vicinity of a hydrated surface the properties of water can differ significantly fromthat of the bulk liquid [1]. These differences are determined by the details of the interfacialenvironment and they are thus sensitive to the specific chemical and topological features ofthe hydrated surface [2]. The details of water’s interfacial molecular structure, therefore,contain information about these surface features and how they affect their local aqueousenvironment [3]. This information is valuable because it provides insight into the collectivemolecular effects that control the solvation properties of complex solutes, but it is also diffi-cult to access due to limitations in our ability to measure and characterize water’s interfacialmolecular structure.Our current understanding of the interfacial molecular structure of liquid water is derivedprimarily from the results of interface sensitive experimental techniques such as vibrationalsum frequency generation spectroscopy [4–8], THz absorption spectroscopy [9, 10], dynamicnuclear polarization [11], NMR [12, 13], and X-ray and neutron scattering [14–16]. Unfor-tunately, these experiments are typically more difficult to interpret than their bulk phasecounterparts due to the constraints associated with achieving interface selectivity. This hasdriven an increased demand for theoretical developments that can facilitate the analysis andinterpretation of these interface-sensitive experiments. Resulting efforts have relied heavilyon the use of atomistic simulation to provide the molecular details of water’s interfacialstructure.Classical molecular dynamics (MD) simulations are a particularly efficient theoreticalframework for modeling the microscopic properties of aqueous interfacial systems. Thesesimulations provide a molecular-level basis for understanding the microscopic structure ofliquid water and how it responds to the anisotropic environment of the liquid phase bound-ary. This response is mediated by the properties of water’s hydrogen bonding network andtherefore involves the correlated arrangements of many individual water molecules. Char-acterizing this high-dimensional molecular structure in simple and intuitive terms can be asignificant challenge, especially for solutes such as proteins that exhibit irregular or hetero-geneous surface properties.Here we address this challenge by characterizing water’s interfacial molecular structurein terms of a low dimensional parameter that quantified its similarity to the structure of2arious interfacial reference systems. We introduce a theoretical framework for quantifyingthis similarity based on the statistical analysis of molecular orientations at the interface. Inthis analysis reference systems serve to designate the unique orientational signatures of waterinterfaces at surfaces with specific well-defined chemical or topographical properties. Theproperties of these reference surfaces can be systematically chosen in order to analyze specificinterfacial features that may be relevant to a particular system of study. This frameworkprovides a physically insightful measure of interfacial structure that eliminates the need toformulate high-dimensional collective variables. Furthermore, by applying this frameworkacross a variety of different reference systems, it can be adapted to report simultaneously onmultiple specific interfacial properties.The general formalism for our framework begins with the definition of the reference sys-tem(s) that will serve as a basis for interfacial characterization. Prior to applying thisframework each reference system must be thoroughly sampled in order to establish its uniqueorientational molecular signature. We quantify this signature in terms of the molecular ori-entational distribution function, f ( (cid:126)κ | ref) = P ( (cid:126)κ | ref) P ( (cid:126)κ | bulk) , (1)where P ( (cid:126)κ | ref) and P ( (cid:126)κ | bulk) denote the equilibrium probabilities for observing a moleculewith a specific molecular configuration, (cid:126)κ , within the given reference system ands in thebulk liquid respectively. We specify molecular configuration in terms of a three-dimensionalvector, (cid:126)κ = (cos θ , cos θ , a ), where θ and θ denote the angles made between the localsurface normal and each of a water molecule’s OH bond vectors, and a denotes the distanceof the water molecule from the instantaneous position of the water interface. We define theinstantaneous water interface following the procedure of Ref. [17]. This system of molecularcoordinates is illustrated schematically in Fig. 1.We characterize the interfacial molecular structure of of a particular interfacial system bysampling values of (cid:126)κ and comparing them to the distribution function, f ( (cid:126)κ | ref). Specifically,we compute the quantity, λ ref = 1 N N (cid:88) i =1 − ln [ f ( (cid:126)κ | ref)] , (2)3 hydrated surface Liquid water ✓ ✓ surface normal FIG. 1. A schematic illustration of the molecular coordinates that are used to specify the orienta-tional state, (cid:126)κ , of molecules at the liquid water interface. where the summation is taken over a set of N molecular configurations sampled from thesystem of interest. This quantity reflects the likelihood for the sampled set of configurationsto occur spontaneously within the environment of the reference system. This likelihood isrelatively large (corresponding to lower values of λ ref ) when interfacial molecular structure issimilar to that of the reference system and relatively small (corresponding to higher valuesof λ ref ) when interfacial molecular structure differs from that of the reference system. Inthis way, λ ref can be used to distinguish hydrated surfaces based on how they influence theirhydration environment, irrespective of their specific chemical or topographical properties.For heterogeneous surfaces, this framework can be applied locally to generate spatiallyresolved maps of λ ref . This can be accomplished by restricting sampling of (cid:126)κ to specificsystematically controllable regions along the surface. A general expression for Eq. (2) thatcan be used to compute λ ref at a specific position, (cid:126)r surf , and time, t is given by, λ ref ( (cid:126)r surf , t ) = 1 τ t + τ (cid:88) t (cid:48) = t − ln [ f ( (cid:126)κ ( (cid:126)r surf , t (cid:48) ) | ref)] , (3)here (cid:126)κ ( (cid:126)r surf , t (cid:48) ) is the orientational configuration of the water molecule that is closest toposition (cid:126)r surf at time t (cid:48) , and the summation is taken over a series of τ discrete time steps. Thevalue of τ in Eq. (3) can be varied to highlight average interfacial response to a heterogeneoussurface ( i.e. , large τ ) or to highlight transient fluctuations in interfacial molecular structure( i.e. , small τ ).Early approaches to mapping the hydration properties of heterogeneous solutes, mostnotably those developed by Kyte and Doolittle, were based on a spatial decomposition4f surface chemistry [18], and have since been extended to provide higher resolution [19].These approaches often fail to accurately predict solvation properties due to their neglect oftransverse correlations within the water interface. More recent approaches have focused onwater-based mapping methods. This includes approaches based on local density fluctuations[20, 21], single water chemical potentials [22], and local electrostatic fields [23]. Our methodis complementary to these previous approaches and can be adapted, via the selection ofdifferent reference systems, to map a wide variety of interfacial properties.We illustrate the performance of our method by applying it to characterize the hydropho-bicity of heterogeneous surfaces. To do this we use a single reference system as a basisfor interfacial characterization – the liquid water interface at an ideal hydrophobic surface.The collective molecular arrangements that are common to this reference system are thusspecified by f ( (cid:126)κ | phob). As we demonstrate below, and in the Supporting Information (SI), λ phob is capable of distinguishing between the interfacial molecular structure of hydrophobicand hydrophilic surfaces and can thus be treated as an order parameter for hydrophobicity.As such, we use λ phob , computed according to Eq. (3), to analyze water’s interfacial re-sponse to heterogeneous surfaces. Since λ phob is based only on water’s interfacial molecularstructure, it can be used to generate hydrophobicity maps that reveal the effective solvationcharacteristics of surfaces with complex or unknown properties.As a proof of concept, we apply our framework for interfacial characterization to a modelsilica surface with a patterned composition of hydrophobic and hydrophilic surface sites[24]. As illustrated in Fig. 2, the surface sites of this model can be either nonpolar ( i.e. ,hydrophobic), if they are terminated with a neutral silica atom, or polar ( i.e. , hydrophilic),if they are terminated with a charged hydroxyl group. Artificial surfaces with well-definedsurface patterns can be created by specifying the hydroxylation state of the surface sites.We then use λ phob to analyze water’s response to various surface patterns. We quantify thisresponse by computing δλ ref = λ ref − (cid:104) λ ref (cid:105) , where (cid:104)· · · (cid:105) denotes an equilibrium averagetaken within the ensemble of configurations sampled directly from the reference system. Inthis way, values of δλ ref ≈ λ phob to distinguish between water’s interfacial molec-5 y SiO O hydrophobic surface site a bc
SiO OO H
Hydrophilic surface site δ + δ − δ + FIG. 2. (a) A snapshot of a simulation of a periodically replicated slab of liquid water in contactwith a 6 × model silica surface. (b) and (c) The chemical termination of the surface sitesdetermines whether they are hydrophobic or hydrophilic. ular response to hydrophobic and hydrophilic regions of a spatially heterogeneous surface.Specifically, we have computed λ phob , using Eq. (3), for points along an extended hydropho-bic surface with a larger rectangular hydrophilic patch, as shown in Fig. 3a. For the dataplotted in Fig. 3b, λ phob has been averaged over a long observation time of τ = 4 ns (40000individual configurations). We observe that water’s average interfacial molecular structureexhibits spatial variations that mimic the chemical patterning of the underlying silica sur-face. Over non-polar regions of the surface δλ phob ≈ i.e. , λ phob ≈ (cid:104) λ phob (cid:105) ), indicatingthat interfacial molecular structure is similar to that of the hydrophobic reference system.Over polar regions of the surface δλ phob > i.e. , λ phob > (cid:104) λ phob (cid:105) ), indicating that water’sinterfacial response differs significantly from that of the reference system. Low amplitudespatial modulations in λ phob can be observed over both the polar and the non-polar regionsof the surface. These modulations reflect the corrugation of the atomic surface and thusindicate that this order parameter is sensitive to the subtle influence of surface topographyon water’s interfacial molecular structure.Transient fluctuations in local interfacial structure can be analyzed by computing δλ phob with a smaller value of τ . For instance, Fig. 3d shows δλ phob computed for the surface in6 bc d x y . . .
52 0 10 20 30 40 50 60 ph o b x/ nm . . . ph o b x/ nm ⌧ = 20ps ⌧ = 2ns y / n m x/ nm . . τ = y / n m x/ nm . . τ = . . . ph o b x/ nm ⌧ = 20ps ⌧ = 2ns . . . ph o b x/ nm ⌧ = 20ps ⌧ = 2ns τ = τ = . . . ph o b x/ nm FIG. 3. (a) A snapshot of the water exposed face of a model silica surface featuring a hydrophilicpatch against a hydrophobic background. (b) and (d) A plot of δλ phob , indicated by shading,computed for points along the surface of the structure shown in Panel (a). The data in Panels (b)and (d) reflects an average over an observation time of τ = 4 ns and τ = 20 ps, respectively. (c) Aplot of the value of δλ phob computed along a line at y = 3 nm (the yellow dotted line in Panel (a))that highlights how interfacial molecular structure is affected by the patch boundary. Grey pointsare values of δλ phob , computed with τ = 20 ps, sampled at different points in time and the solidline is δλ phob computed with τ = 4 ns. Fig. 3a using a value of τ = 20 ps (100 individual configurations). The use of a shorterobservation time highlights the transient details of water’s interfacial molecular structure.Thus, by comparing local values of δλ phob over multiple consecutive snapshots it is possibleto observe the transient fluctuations of interfacial molecular structure and investigate howthey depend on the details of local surface chemistry and topology. The dynamic range of δλ phob within the hydrophobic reference system is described in the SI.Water molecules that reside over the boundaries between polar and non-polar regionsof the surface experience a laterally anisotropic aqueous environment. In these regions,such as along the edge of the polar patch of the surface illustrated in Fig. 3a, δλ phob δλ phob in these boundary regions reveal details about themolecular correlations that mediate interactions along liquid water interfaces, and how thesecorrelations are influenced by the details of surface-water interactions.In the case of the model silica surface we observe that the effect of the polar/non-polarsurface boundary on the interfacial molecular structure is local, limited to the region directlyabove the surface boundary. A cross-section of δλ phob that cuts through the center of thehydrophilic surface patch is plotted in Fig. 3c. This plot reveals that the influence of alarge hydrophilic surface patch on water’s interfacial molecular structure only extends aboutone molecular diameter beyond the edge of the patch. Evidently, in this case any long-ranged transverse correlations due to heterogeneous surface chemistry are not supported bydistortions of the interfacial molecular structure. We present results for a variety of differentsurface patterns in the supporting information (SI).The results presented in Fig. 3 verify that δλ phob is effective as a local order parameterfor surface hydrophobicity. This order parameter can thus be used to infer the effectivehydration properties of an unknown aqueous surface. For the model silica surfaces theheterogeneity in surface chemistry is mirrored by the spatial dependence of δλ phob , however,for more complex surfaces, the relationship between surface structure and δλ phob is not sostraightforward. When this is the case, δλ phob can provide valuable physical insight into therelationship between surface heterogeneity and local hydration properties.The complex heterogeneous surface properties of hydrated proteins are reflected in thewater’s spatial dependence on interfacial molecular structure. Figure 4 illustrates that thisspatial dependence can be revealed with δλ phob . Specifically, Fig. 4 illustrates the valueof δλ phob computed along the surface of the inactive CheY protein (PDB code: 1JBE) [25]using Eq. (3) and a value of τ = 10 ps. This map of δλ phob indicates regions of the proteinsurface whose interactions with water result in hydrophobic ( i.e. , green shaded regions) orhydrophilic ( i.e. , purple shaded regions) interfacial molecular structure.Details of these protein surface maps, such as the position and shapes of the hydrophilicdomains, are sensitive to the conformational fluctuations of the protein. This is illustratedin Fig. 4c, which shows the map of δλ phob computed for the CheY protein for a consecutive8 bc -0.5100.5 δλ phob s u r f a ce r e s i du e t/ ns . . . . . . Time
FIG. 4. (a) A simulation snapshot of the CheY protein. Water molecules are omitted for clarity. (b)A map of δλ phob computed for points along the surface of the protein using τ = 10 ps. Color scaleis identical to that in Fig. 3b,d. (c) Spatial maps of δλ phob for a series of protein conformationsspaced out along a 4 ns trajectory. sequence of conformations. The dynamics of interfacial structure can be further analyzedby considering the time dependence of δλ phob for individual surface residues. Figure 5 showsthe result of such a calculation performed along a 10 ns trajectory of the CheY protein.We observe that the interfacial structure in some regions of the protein remains relativelystatic, as indicated by persistent green or purple bands in Fig. 5a, while other regions ofthe protein exhibit interfacial structure that fluctuates significantly in response to proteinconformational dynamics.The ability to map the hydration dynamics of fluctuating irregular solutes is a uniquefeature of this interfacial characterization method. Of course, this method is not limitedto the applications presented above. For instance, the approach can be easily extended toreport on additional hydration properties with the use of different reference systems, such assystematically charged surfaces or those with specific curvature, or by expanding the defini-tion of (cid:126)κ , for example to include dynamical information. Nor is the general approach limitedto extended liquid water interfaces. The examples described here simply demonstrate thetype of insight that can be derived from analyzing the orientational properties of interfacialwater molecules. 9 s u r f a ce r e s i du e t/ ns . . . . . . s u rf ace r e s i du e s u r f a ce r e s i du e t/ ns . . . . . . . . . . . . . . . ph o b t/ ns -0.200.20.40.6 δ λ phob Res δλ phob ab FIG. 5. (a) A time series plot of δλ phob , indicated by shading, as computed for individual surfaceresidues of the CheY protein. In this plot, each row corresponds to a unique surface residue. (b)A plot highlighting of the dynamics of δλ phob for three specific surface residues. ACKNOWLEDGEMENT
We acknowledge useful discussions with Professor John Weeks and Professor GerhardHummer. This work was supported by the National Science Foundation under CHE-1654415and also partially (SS) by the Kwanjeong Educational Foundation in Korea. ∗ [email protected][1] K. B. Eisenthal, Chemical reviews , 1343 (1996).[2] D. Qu´er´e, Annual Review of Materials Research , 71 (2008).[3] A. Verdaguer, G. M. Sacha, H. Bluhm, and M. Salmeron, Chemical reviews , 1478 (2006).[4] Q. Du, R. Superfine, E. Freysz, and Y. R. Shen, Physical Review Letters , 2313 (1993).
5] M. Schleeger, Y. Nagata, and M. Bonn, The Journal of Physical Chemistry Letters , 3737(2014).[6] P. C. Singh, The Journal of chemical physics , 094706 (2012).[7] S. Baldelli, C. Schnitzer, and D. Simonelli, The Journal of Physical Chemistry B , 5313(2002).[8] D. E. Gragson, B. M. McCarty, and G. L. Richmond, Journal of the American ChemicalSociety , 6144 (1997).[9] V. Conti Nibali and M. Havenith, Journal of the American Chemical Society , 12800(2014).[10] S. Ebbinghaus, S. J. Kim, M. Heyden, X. Yu, U. Heugen, M. Gruebele, D. M. Leitner, andM. Havenith, Proceedings of the National Academy of Sciences , 20749 (2007).[11] B. D. Armstrong and S. Han, Journal of the American Chemical Society , 11270 (2009).[12] G. Otting and K. Wuethrich, Journal of the American Chemical Society , 1871 (1989).[13] G. Otting, E. Liepinsh, and K. W¨uthrich, Science (New York, N.Y.) , 974 (1991).[14] D. I. Svergun, S. Richard, M. H. J. Koch, Z. Sayers, S. Kuprin, and G. Zaccai, Proceedingsof the National Academy of Sciences , 2267 (1998).[15] F. Bruni, M. A. Ricci, and A. K. Soper, The Journal of chemical physics , 1478 (1998).[16] G. Luo, S. Malkova, J. Yoon, D. G. Schultz, B. Lin, M. Meron, I. Benjamin, P. Vanysek, andM. L. Schlossman, Science (New York, N.Y.) , 216 (2006).[17] A. P. Willard and D. Chandler, The Journal of Physical Chemistry B , 1954 (2010).[18] J. Kyte and R. F. Doolittle, Journal of molecular biology , 105 (1982).[19] L. H. Kapcha and P. J. Rossky, Journal of molecular biology , 484 (2014).[20] D. T. Limmer and A. P. Willard, Chemical Physics Letters , 144 (2015).[21] A. J. Patel, P. Varilly, D. Chandler, and S. Garde, Journal of statistical physics , 265(2011).[22] T. Young, R. Abel, B. Kim, B. J. Berne, and R. A. Friesner, Proceedings of the NationalAcademy of Sciences , 808 (2007).[23] R. C. Remsing and J. D. Weeks, The Journal of Physical Chemistry B , 9268 (2015).