Computation of Overhauser dynamic nuclear polarization processes reveals fundamental correlation between water dynamics, structure, and solvent restructuring entropy

Abstract

Hydration water dynamics, structure, and thermodynamics are crucially important to understand and predict water-mediated properties at molecular interfaces. Yet experimentally and directly quantifying water behavior locally near interfaces at the sub-nanometer scale is challenging, especially at interfaces submerged in biological solutions. Overhauser dynamic nuclear polarization (ODNP) experiments measure equilibrium hydration water dynamics within 8–15 angstroms of a nitroxide spin probe on instantaneous timescales (10 picoseconds to nanoseconds), making ODNP a powerful tool for probing local water dynamics in the vicinity of the spin probe. As with other spectroscopic techniques, concurrent computational analysis is necessary to gain access to detailed molecular level information about the dynamic, structural, and thermodynamic properties of water from experimental ODNP data. We chose a model system that can systematically tune the dynamics of water, a water–glycerol mixture with compositions ranging from 0 to 0.3 mole fraction glycerol. We demonstrate the ability of molecular dynamics (MD) simulations to compute ODNP spectroscopic quantities, and show that translational, rotational, and hydrogen bonding dynamics of hydration water align strongly with spectroscopic ODNP parameters. Moreover, MD simulations show tight correlations between the dynamic properties of water that ODNP captures and the structural and thermodynamic behavior of water. Hence, experimental ODNP readouts of varying water dynamics suggest changes in local structural and thermodynamic hydration water properties.

---

Introduction

Local hydropathies near molecular interfaces to water (e.g., proteins and polymers) modulate the surface activity for solute binding,^1–5 among many other properties.^2,6–19 Locally hydrophobic regions of proteins are key to facilitating folding and inter-protein interactions, and are often characteristic of active sites.^1,16,17,20 In reality, hydropathy or hydrophobicity is not truly a property of the protein itself, but rather a statement of how water in the vicinity responds to the presence of particular protein side chains. Experimentally characterizing the structural, dynamical, and thermodynamic properties of water in fully hydrated environments at the molecular scale is challenging. Several experimental methods can probe the behavior of hydration layer waters,^{6,16,19,21,22} but experimental measurements are often limited to detecting average water properties of the entire ensemble.^23–25 There is extensive literature on molecular dynamics (MD) simulation studies that probe hydration dynamics of different solution systems with heterogeneous water–protein interfaces,^{7–9,11,16,17,21,26–34} but relatively few directly compare and validate with experimental techniques,^16,17,21,22 in part, due to the difficulty of resolving local water dynamics experimentally. Undoubtedly, it is critical to leverage atomistic MD simulations synergistically with experiments to comprehensively and robustly characterize the heterogeneous hydration environments on soft material and biomolecular surfaces^{1,3,4,9,17,20}

Advanced techniques such as quasi-elastic neutron scattering (QENS),^17,35,36 Overhauser dynamic nuclear polarization (ODNP),^16,21,37 terahertz (THz) spectroscopy^38–41 and pump–probe infrared (IR) spectroscopy⁴² can directly yield surface-specific water properties. Among those, only ODNP is sensitive to translational hydration water properties around localized sites or surfaces that are fully surrounded by bulk water in biological or synthetic solution-state environments. Fundamentally, ODNP is a magnetic resonance technique that quantifies electron-¹H cross-relaxation by measuring (1) the enhancement of ¹H NMR signals induced by dynamic nuclear polarization (DNP) via the transfer of polarization from an unpaired electron of a nitroxide spin probe to the ¹H nucleus of water, and (2) the longitudinal spin lattice relaxation time, T₁, that reports on all ¹H relaxation mechanisms induced by the dipolar coupling between the ¹H nuclear spin and the electron spin of the nearby spin probe. The electron-¹H cross-relaxivity at a magnetic field of 0.35 Tesla and electron Larmor frequency of 9.8 GHz is sensitive to translational movement of hydration waters near the electron spin probe on instantaneous timescales ranging from about 10 picoseconds to nanoseconds and within 8 to 15 angstroms of nitroxide-based spin probes that can be tethered to specific surface sites.

ODNP responds sensitively to the chemically or topologically distinct local environments that are generated by highly heterogeneous surfaces immersed in bulk water. It thus demonstrates a diversity of local water dynamics and structure that reflects the diversity in these local environments. A recent synergistic ODNP-MD simulation study of sites with varying hydropathy on a CheY protein surface displayed a positive correlation between site-specific hydrophobicity and translational water dynamics.^16,43 This study demonstrated a connection between the ODNP spectroscopic quantities and computed thermodynamic properties of water by demonstrating a correlation between the two. However, no study has directedly computed ODNP parameters for the purpose of deriving hydration water dynamical information and coupling to molecular dynamical, structural, and thermodynamic properties of hydration water near surfaces or solutes.

In the present work, we compute ODNP spectroscopic parameters—specifically, the coupling factor and electron-¹H cross-relaxation rates—and use molecular dynamics simulations to characterize the dynamical, structural, and thermodynamic properties of water in water–glycerol mixtures with increasing glycerol content [Fig. 1]. Water–glycerol mixtures allow us to probe the complex interrelationship between water properties at the molecular scale while varying solution viscosity by a known amount. Moreover, glycerol is of fundamental interest due to its ubiquity in biological studies for its role in cryopreservation of proteins and has been recently shown to alter not only the dynamical, but also structural and thermodynamic properties of water.⁴⁴ Regardless of the system, the study of molecular determinants of surface hydration requires dual experimental and computational insight. Experimental techniques like ODNP offer valuable, but incomplete, insight on the molecular scale properties of hydration waters. On the other hand, without experimental validation, simulations can be subject to doubts, especially when it comes to variations in the structure and dynamics of water near interfaces. Synergistic fully-atomistic simulations can aide in elucidating the molecular details that are often inaccessible to experimental methods.

Fig. 1 Snapshots of simulation boxes for a range of glycerol–water compositions. Here, OPC water molecules are represented as VDW spheres and glycerol is shown in green (licorice representation). 4-OH-TEMPO is shown in licorice representation with carbons, oxygens, hydrogens, and nitrogens represented in cyan, red, white, and blue, respectively. The water and glycerol molecules closest to the 4-OH-TEMPO are made fully transparent in each panel to make the 4-OH-TEMPO visible amidst the densely packed glycerol mixture.

In this study, we perform atomistic MD simulations of glycerol–water to test the ability of MD simulations to reproduce ODNP spectroscopic quantities and to connect ODNP spectroscopic and MD-derived translational, rotational, and hydrogen bonding dynamical quantities. We furthermore exploit the atomistic information content of MD simulations to quantify the relationship between hydration water dynamics, structure, and thermodynamics in glycerol–water mixtures across a wide range of compositions.

Methods

Overhauser dynamic nuclear polarization

We perform experimental ODNP measurements with samples containing mixtures of water (Milli-Q purity, Milli-Q UV Plus system, Millipore Inc., Bedford, MA, USA), glycerol (d₈, 99%, DLM-558, Cambridge Isotope Laboratory, Tewksbury, MA, USA) and 4-hydroxy-TEMPO (176141, Sigma Aldrich) at a constant concentration of 318 μM. We transfer 4 μL samples into round quartz capillaries of 0.6 mm ID × 0.84 mm OD (Vitrocom, New Jersey, USA), sealed with Critoseal on one end and melted beeswax on the other end. We perform ODNP experiments using the X-channel of a a Bruker EMXPlus spectrometer and a Bruker Avance III NMR console (Bruker, Massachusetts, USA) that utilizes specialized automation (AU) programs to interact with the ODNP microwave hardware. We install the capillaries in a home-built NMR probe with a U-shaped coil centered in a high sensitivity microwave cavity (Bruker ER 4119HS-LC).⁴⁵ ODNP measurements rely on the saturation of the electron paramagnetic resonance (EPR) of 4-hydroxy-TEMPO at 9.8 GHz (X-band) frequency at a center field of 3484 G, which coincides with the center resonance of the nitroxide EPR spectrum. Dry air at 9.44 L min⁻¹ streams through the cavity across the probe and capillary to maintain a temperature of 18 °C. Theory and experimental details are outlined in previous publications.^16,22,37,46

Molecular dynamics simulations

We perform molecular dynamics simulations using the OPC 4-site water model,⁴⁷ the Blieck–Chelli (BC) model for glycerol^48–50 and a 4-hydroxy-TEMPO nitroxide spin probe. Both the OPC water model and the BC model for glycerol accurately reproduce the diffusivities of pure water⁴⁷ and pure glycerol^48,50 under ambient conditions (298.15 K and 1-bar). For the spin probe, its partial charges are obtained using the AMBER18 Antechamber package⁵¹ informed by quantum chemical calculations using the Gaussian 16 software.⁵² Specifically, we apply the B3LYP functional and the 6-311++G(d,p) to perform geometry optimization of 4-hydroxy-TEMPO. All other inter- and intramolecular parameters derive from the second-generation generalized Amber forcefield (GAFF2).^53,54 The results of this parametrization scheme yield similar parameters to those obtained in previous publications.⁵⁵ All Coulombic interactions are modeled with the particle–mesh Ewald summation scheme (PME).⁵⁶

We simulate glycerol–water–spin probe systems with glycerol mole fractions (x_glyc) of 0, 0.01, 0.033, 0.05, 0.075, 0.1, 0.15, 0.2 and 0.3 using the GPU-optimized OpenMM molecular simulation software.⁵⁷ At each of these concentrations, we include a single 4-OH-TEMPO molecule in the simulation box. Hence, the spin probe concentration varies depending on glycerol concentration (from 27 mM in pure water to 10 mM for x_glyc = 0.3). We first energy minimize each system, then equilibrate in the NPT ensemble using a Langevin thermostat⁵⁷ paired with a Monte Carlo barostat⁵⁷ at 290 K and 1 atm. Following equilibration, the NPT run continues for 250 ns with system configurations saved every 10 ns. Each saved configuration serves as the starting point for an independent 1 nanosecond NVE simulation for dynamic properties, with system coordinates saved every 0.1 ps. To calculate hydration water dynamics timescales and ODNP spectroscopic quantities, we compute 95% confidence intervals by bootstrapping the results obtained from 20 independent MD simulations.

We characterize the effect of glycerol on solvation thermodynamics via the solvation free energy ΔG_solv of methane—an ideal small hydrophobic solute—in glycerol–water. To estimate ΔG_solv for this series of mixtures, we implement an expanded ensemble simulation procedure in which we gradually scale intermolecular interaction parameters between the methane molecule and the glycerol–water mixture. We smoothly scale Lennard Jones (LJ) and coulombic interaction parameters via a scalar parameter λ from λ = 0 (non-interacting, or ideal gas molecule) [panels (1) and (3) in Fig. 6(a)] to 1 (fully interacting methane) [panels (2) and (4) in Fig. 6(a)]. To estimate the Gibbs free energy of solvation ΔG_solv [Fig. 6(b)], we then apply the multistate Bennett acceptance ratio (MBAR) method distributed via the pymbar⁵⁸ Python library.

Theory

Computing spectroscopic quantities

Dipolar autocorrelation functions. ODNP is an NMR technique that quantifies pairwise magnetic dipolar cross-relaxation between an electron spin on a free radical spin probe molecule (here of 4-hydroxy-TEMPO) and the nuclear spin of the water proton. The spin–spin dipolar coupling energy between the radical electron and a water proton that are dynamically diffusing in solution state is governed by the following semi-classical Hamiltonian^24,59where F^(m)₂(r(t)) (and F^(m)₂(r(t))*) are spherical harmonic functions (and their complex conjugates) [ESI,† Section D] that are dependent on the displacement vector between the electron spin and the proton nuclear spin, r(t). Î_i and Ŝ_i are the quantum mechanical spin operators for the proton nuclear spin and electron spin, respectively. The constant c_δ = μ₀ħγ_Sγ_I/4π appears in eqn (1) where , , and γ_S = γ_I/1.52 × 10⁻³ are the vacuum permeability, gyromagnetic ratio of the proton spin, and the gyromagnetic ratio of the electron spin, respectively.

Using this semi-classical framework, we compute the equilibrium translational diffusion of water molecules by constructing the time-autocorrelation functions (ACFs), C^(m)_ODNP(t), that depend purely on the relative positions of water hydrogen and the unpaired electron of a spin probe via the classical spherical harmonic functions included in eqn (1)

where is the displacement vector between the water hydrogen i and the oxygen radical of the spin probe and m = 0, 1, 2 the order, which correspond, respectively, to the flip–flop, single spin flip, and double spin flip transitions of the coupled I–S spin system. The C^(m)_ODNP(t) functions decay with time as water molecules diffuse from their initial position relative to the spin probe, just as the dipolar coupling between the proton and electron spin weakens with increased spatial separation. The C^(m)_ODNP(t) functions are complex-valued, but the complex part is negligible for isotropic systems⁵⁹ and thus ignored in this study [more detailed description in ESI,† Section D]. Further, the system isotropy ensures that C^(m)_ODNP(t) = C⁽ⁿ⁾_ODNP(t) for all n and m. Hence, we refer to only a single ACF for computing ODNP properties, C_ODNP(t), for the remainder of the discussion. For the present work, we find that all measured C_ODNP(t) are well-described by a tri-exponential fit:

C_ODNP,fit(t) = a₁e^−t/τ₁ + a₂e^−t/τ₂ + a₃e^−t/τ₃ (3)

where τ₁ > τ₂ > τ₃ and the coefficients, a_i (i = 1, 2, 3), sum to unity. Note, that the fitting parameters are independently determined for each glycerol–water mixture.

Spectral density function. With analytical models for C_ODNP(t) at hand [eqn (3)], we next derive the spectral density function, J(ω), required to compute ODNP spectroscopic quantities for a given glycerol–water mixture. Specifically, we define the three spectral densities via the real part of the Fourier transform of C_ODNP,fit(t) that contains three decay terms, and hence three τ values.Here, the three terms in the sum account for the contribution of long (τ₁), intermediate (τ₂) and short (τ₃) timescales to local water diffusion. In the low-frequency regime (ω ≪ ω_S), the spectral densities are more affected by the long timescale contribution (collective motion). On the other hand, the short timescale contribution (instantaneous motion) dominates the spectral densities at high frequency (ω ≫ ω_I). To computationally derive the spectroscopic quantities measured by ODNP, we determine the transition rate via the amplitude of the spectral density, J(ω), for an electron spin–proton spin (e–n) transition that occurs at a given electron spin resonance (ESR) frequency, ω. Alternatively, J(ω) may be predicted using a theoretical continuum model⁶⁰—the force-free hard sphere (FFHS) model—for water diffusion relative to the spin probe's radical electron. While previous work has demonstrated that the FFHS approach yields experimental water dynamics in qualitative agreement with simulation results,^21,61 the approach presented in eqn (4) more directly determines the spectral density from e–n dipolar correlations with molecular detail and in dynamic systems that cannot be modeled by the FFHS model.

Fig. 2(a) illustrates the dynamic spin–spin dipolar interaction between the unpaired electron of 4-OH-TEMPO and water protons underpinning the calculation of C_ODNP,fit(t). In Fig. 2(b), we depict C_ODNP(t) for a representative model system, a water–glycerol mixture with x_glyc = 0.1. We find that C_ODNP(t) and its model fit, C_ODNP,fit(t), are nearly indistinguishable with R² = 0.99. In Table S1 (ESI†), we summarize the tri-exponential fitting parameters for each glycerol–water mixture. In Fig. 2(c), we apply eqn (4) to compute the spectral density function J(ω) from C_ODNP(t) for the same mixture with x_glyc = 0.1. Here, the amplitude of the spectral density values at frequencies ω_I and ω_S [J(ω_I = 14.8 MHz) and J(ω_S = 9.8 GHz)] are sensitive to translational dynamics on nanosecond and picosecond timescales, respectively.

Fig. 2 Schematic of ODNP spectroscopic quantities calculation from classical MD trajectories. (a) The snapshot shows 4-OH-TEMPO in a 0.1 mole fraction glycerol mixture at 290 K with bulk water represented in VMD as a medium and the hydration waters within 3.5 angstroms of the spin probe in VDW sphere representation. Glycerol molecules are omitted for clarity. The diagram to the right of the snapshot illustrates the ODNP mechanism for the nearest water molecule to the spin probe at some time t. (b) The ODNP time autocorrelation function C_ODNP(t) at a glycerol mole fraction of 0.1 (red line) is fit to a tri-exponential model [eqn (3)] (black line) as described in the text. (c) The real part of the Fourier transforms of C_ODNP,fit(t) gives spectral density function J(ω). Subsequently, J(ω) values at radical electron J(ω_S) and proton Larmor frequencies J(ω_I) are identified on the plot by the blue and green vertical lines, respectively. Approximate functional forms of ODNP spectroscopic quantities k_σ and ξ are in blue and green text, respectively. (d) Based on the force-free hard sphere model,⁶⁰ we depict the expected trends in k_σ/k_σ,bulk (blue line) and ξ/ξ_bulk (green line). Here, k_σ,bulk and ξ_bulk are the expected bulk values of k_σ and ξ in pure STP water as defined by Franck et al.²²

ODNP spectroscopic probes. From the amplitudes of the spectral density function, we directly compute several spectroscopic quantities derived from ODNP experiments. One such quantity is the cross-relaxivity k_σ—the rate of transitions for the mutual flip of electron and proton spins in opposite directions—that strongly depends on the translational mobility of water on instantaneous timescales and within 8–15 angstroms of the spin probe.^16,37,60 As defined in prior work,^46,60,62k_σ originates from the zero and double quantum transition of the dipolar coupled electron and proton spin pairs (expressed in the Hamiltonian of eqn (1)), as given bywhere C_SL, ω_S = γ_SB₀, and ω_I = γ_IB₀ are the molar spin label concentration [Table S2, ESI†], Larmor precession frequency of the electron spin, and Larmor precession frequency of the proton spin, respectively. The ODNP experiments are conducted at a static magnetic field of B₀ = 0.35 T, which sets ω_S = 9.8 GHz and ω_I = 14.8 MHz. For the purposes of the present study, we determine the k_σ values from the simulation-derived spectral densities [eqn (4)].

To provide an illustration, we depict the theoretical FFHS-derived cross-relaxivity, k^FFHS_σ, at B₀ = 0.35 T as a function of the correlation time in Fig. 2(d). The k^FFHS_σ values exhibit non-monotonic behavior with increasing translational correlation time, τ_c, of hydration waters. Specifically, k^FFHS_σ increases with correlation times up to τ_c ≈ 100 ps. Upon further retardation of translational dynamics, k^FFHS_σ decreases monotonically with increasing τ_c. Notably, this maximum in k^FFHS_σ is found at τ_c that is 3 times greater than the expected correlation time, τ_c ≈ 33 ps, for pure water at standard temperature (298.15 K) and pressure (1 bar) conditions (STP).²² The precise location of the peak in k^FFHS_σ strongly depends on the Larmor precession frequencies by way of the applied magnetic field strength, B₀.

Under the present B₀, hydration waters with translational diffusion coefficients smaller than 1/3 of the diffusivity of pure water under ambient conditions (T = 25 °C and P = 1 bar), , yield k_σ values that linearly decrease with the diffusivity of hydration waters, D_H2O,local. In systems with , such as waters hydrating catalyst support surfaces, hydrophilic materials interfaces,⁴² or protein surfaces and interiors,¹⁶k_σ decreases upon retardation of hydration water translational dynamics. However, under the conditions considered in the present work, k_σ begins to decrease when . Here, we demonstrate the non-monotonicity of k_σ by spanning water translational dynamics from those of pure water to those of water–glycerol mixtures with .

Another ODNP spectroscopic quantity of interest is the proton self-relaxivity k_ρ—the rate at which water proton polarization returns to thermal equilibrium modulated by the electron spin dipolar coupled to the proton nuclear spin, as follows:

Unlike k_σ, k_ρ is the rate of transition for all proton nuclear spin flip events induced by dipolar coupling to the electron spins, not only the mutual proton–nuclear spin flips, and hence depends on the spectral density at the sum and difference of the electron and proton Larmor frequency and the proton Larmor frequency alone. The amplitude of the spectral density at the sum and difference frequencies, ω_S ± ω_I, is close to simply the electron spin Larmor frequency, ω_S, and is sensitive to instantaneous translational dynamics correlation times of sub-nanoseconds, while the amplitude of the spectral density at the nuclear Larmor frequency, ω_I, is sensitive to longer 1–10 ns timescales for the slower, collective, motion of water. Because both k_σ and k_ρ contain the pre-factor the coupling factor that is defined by the ratio of k_σ and k_ρ eliminates the dependence on C_SL.In experimental systems with difficult-to-quantify spin label concentration, the elimination of the pre-factor in eqn (7) facilitates direct comparison between the experimental and MD-derived ξ. Conveniently, ξ is readily obtained from ODNP measurements of ¹H NMR signal enhancements and ¹H T₁ spin lattice relaxation times.^37,63,64 As with k_ρ, ξ depends on both instantaneous and collective hydration water dynamics.²² When assuming simple diffusion as reflected in the FFHS model to determine the form of the spectral density function,^22,60 the value of ξ^FFHS monotonically depends on the characteristic correlation time for the translational diffusion dynamics of hydration water, simplifying its analysis, and as depicted in Fig. 2(d).

Proton spin–lattice relaxation times. We computationally extract the water proton spin–lattice relaxation times T₁₀[0], which depend on both the system-average translational and rotational dynamics of water. We directly compute T₁₀[0] according to the relationship presented by Bloembergen, Purcell and Pound⁶⁵ which assumes that effects of J-coupling, spin-rotation, and chemical shift anisotropy of the ¹H NMR signals are negligible, and that T₁ only depends on the system-average dynamics of water that is strongly modulated by the viscosity in a glycerol–water mixture.⁶⁵ Because the primary relaxation mechanism involved in T₁₀[0] is dipolar coupling between water protons modulated by water dynamics, we again compute the time autocorrelation functions, this time between the water protons . The right-hand-side of the equation is nearly identical to eqn (2), but with the proton–proton displacement vector replacing . To efficiently compute , we only consider displacement vectors between a randomly chosen “probe” water proton and all other nearby water protons that the probe water proton encounters. Just as with C_ODNP(t), the decays monotonically as water protons diffuse away from the “probe” proton. As in eqn (4), we fit eqn (6) to a multiexponential model and analytically Fourier transform the multiexponential functions to obtain the spectral density functions, as follows:where b^(m)₁, τ₁ and τ₂ (with τ₁ > τ₂) are fit parameters and b₂^(m) = 1 − b^(m)₁. We apply a bi-exponential model here because decays much faster than C_ODNP(t), allowing us to fit with R² > 0.95 for the entire range of glycerol concentrations with only two exponentials, as shown in Table S2 (ESI†). Given these analytical spectral densities, we directly compute the longitudinal relaxation time using a formula from Bloembergen and coworkers⁶⁵Further, we physically interpret the translational and rotational contributions to T₁₀[0] via the decomposition T₁₀⁻¹[0] = (T₁₀⁻¹)_inter + (T₁₀⁻¹)_intra. Here, (T₁₀⁻¹[0])_inter and (T₁₀⁻¹[0])_intra refer to the intermolecular and intramolecular contributions to T₁₀⁻¹[0], respectively. The intramolecular contribution to the relaxation rate, (T₁₀⁻¹[0])_intra, can be ignored when using a rigid water model (such as the OPC 4-site model used in this study) that depends only on the system-average rotational dynamics of water given that the distance between the two hydrogen atoms is fixed. For the intermolecular contribution to (T₁₀⁻¹[0])_inter between water, the spherical harmonic functions strictly depend on the displacement between water protons on separate water molecules, and not the re-orientation of the displacement vector between protons on the same water molecules. As such, reduces to the same functional form as C_ODNP(t) with water protons acting as probes (such as the electron in a nitroxide spin probe) for the surrounding water protons. Because water protons are indistinguishable from each other, the resulting correlation functions are manifestly system-average properties of water. The corresponding intermolecular relaxation rate, (T₁₀⁻¹[0])_inter, calculated from eqn (8) and (9) is hence sensitive to the system-average translational dynamics of water.

Computational probes of water dynamics

Translational dynamics. We characterize the translational dynamics of water by computing the so-called survival probability C_survival(t)^33,61 ACF that quantifies the timescale for waters to remain near the spin probe:where N_w gives the number of water molecules in the simulation box, and S_i(t) an indicator function that is 1 if the molecule i is inside a cutoff radius of 8-Å from the unpaired electron of the spin probe—approximately the width of the first two hydration shells around the radical oxygen. We apply absorbing boundary conditions such that only the water molecules that remain continuously within the cutoff radius from the initial time t = 0 to time t contribute to C_survival(t). We find that neither choosing smaller cutoff radii (for instance, the first hydration shell near 5-Å) nor removing the absorbing boundary conditions qualitatively affect the trends in any ACF described here [Fig. S8, ESI†].

Rotational dynamics. To quantify the rotational dynamics of hydration waters, we compute the orientational ACF (OACFs)^{27–29,31,66} that measures a characteristic time for water reorientation:where N_w gives number of waters within the 8-Å cutoff radius at initial time t = 0, P_l(·) is the l-th Legendre polynomial function, and is the unit vector of water dipole i. For this present work, we only consider the second order OACF C⁽²⁾_OACF due to its relevance to longitudinal spin-relaxation rates.

Hydrogen bond dynamics. To probe the dynamics of water–water hydrogen bonding, we compute the hydrogen bond survival probability³⁰where N_HB is the number of water–water hydrogen bonds containing waters within the cutoff radius at initial time t = 0, and h_i(t) is a function that assumes a value of 1 if hydrogen bond i is intact at time t. We define hydrogen bonds via the widely-used geometric criteria of Luzar and Chandler,⁶⁷ namely, distance and angular cutoff values of 3.5 angstroms and 120 degrees, respectively [see the inset schematic in Fig. 4(c)].

Estimating relaxation time constants. To quantify the shifts in water's equilibrium dynamics with varying x_glyc, we compute several relaxation time constants: the ODNP derived translational diffusion correlation time (τ_ODNP), survival correlation time due to translational diffusion (τ_survival), rotational diffusion correlation time (τ_OACF), and hydrogen bond correlation time (τ_HB). More specifically, we estimate these time constants by integrating the multiexponential fits to all ACFs detailed abovewhere i = ODNP, survival, OACF, or HB.

Results

Direct comparison of ODNP and MD-derived spectroscopic quantities

To directly probe the effect of hydration water retardation on the ODNP parameters (including T₁₀[0]), we systematically increase the solution viscosity by adding glycerol to water at mole percentages ranging from x_glyc = 0 to 0.3. In Fig. 3(a), we highlight the composition-dependent amplitude of the spectral density, J(ω), at the Larmor precession frequencies of protons (ω_I = γ_IB₀ ≈ 14.8 MHz) and unpaired electrons (ω_S = γ_SB₀ ≈ 9.8 GHz) at B₀ = 0.35 T. The increase in J(ω_I) from pure water to x_glyc = 0.3 reflects on the increase in the relaxation rate of the proton spins, which according to eqn (7), leads to a monotonic decay of the coupling factor, as depicted in Fig. 3(b). The inset of Fig. 3(a) shows the approximate spectral density contribution to the cross-relaxation rate of 5J(ω_S). The approximation stems from the limit of ω_S ≫ ω_I and hence—by eqn (5)—k_σ goes as . 5J(ω_S) shows the same trend as the direct computation of k_σ [Fig. 3(c)], increasing up to a glycerol concentration of x_glyc = 0.1 and decreasing at x_glyc > 0.1.

Fig. 3 ODNP spectroscopic quantities measured experimentally and computed from classical MD simulations. (a) ODNP spectral density functions as a function of glycerol concentration where lighter-colored lines correspond to higher glycerol content with the Larmor precession frequency of the proton and radical electron indicated by the green and blue vertical lines, respectively. Comparing experimentally and computationally determined (b) coupling factor ξ, (c) cross-relaxivity k_σ, and (d) T₁₀[0] as a function of glycerol content. ODNP experiments and MD simulations yield spectroscopic quantities with similar trends with increasing glycerol concentration.

The MD simulation-derived relative coupling factor ξ_r = ξ/ξ_pure exhibits a similar monotonic decrease as the experimental coupling factor with increasing glycerol concentration throughout the whole range of the glycerol–water mixture [Fig. 3(b)]. Here, we normalize ξ by the coupling factor of pure water to better compare directly to the experimental results. Though simulation systematically underestimates the experimental ξ_r values by 5 to 60%, we observe qualitative agreement between both measures of ξ_r for x_glyc < 0.30. Notably, the 95% confidence interval (CI) broadens with increasing glycerol concentration. At x_glyc = 0.3, MD simulations underestimate ξ_r by 60%, with the experimental value lying outside the 95% confidence interval (CI) of the MD-computed value. While ODNP measurements measure e–n dipolar correlations averaged over 10¹⁴ spin probes, our simulated systems contain a single spin probe molecule. Hence, we observe inherent fluctuations in C_ODNP(t) due to limited sampling of the dipolar interactions between the spin probe and the water nuclei. Increasing the concentration of glycerol exacerbates this sampling limitation due a systematic decrease in the average number of hydration waters. Increasing the NVE simulation length computes C_ODNP(t) with improved resolution, which we assess by performing 20 additional independent 3-ns long (versus the original 1-ns) simulations at x_glyc = 0.30. For x_glyc > 0.3, the difficulty of sampling water–TEMPO interactions combined with slowing translational dynamics makes accurate calculation of ODNP quantities intractable. While increasing NVE simulation time slightly decreases the width of the 95% CI on ξ_r (approximately by 6%), the simulations continue to systematically underestimate ξ_r at x_glyc = 0.3. This suggests that the disagreement between experiment and simulation ξ_r does not solely stem from under sampling of long timescale collective water motions.

Instead, the quantitative disagreement between the experimental and simulation results may stem in large part from the use of classical, fixed-charge molecular models. For instance, such MD models can only approximate interactions between water and a fixed position atomic position of 4-OH-TEMPO. In reality, the electron spin is delocalized between the N and O of the nitroxide radical. Additionally, the water hydrogen may polarize differently near the spine probe compared to bulk water, yielding different O–H bond length and H–O–H angle. Further, OPC is a rigid water model and hence does not contain bonded (O–H) or angular (∠H–O–H) interaction terms. In particular, the inability to capture O–H vibrations with rigid water molecules likely affects the measurement of proton–nitroxide pair distances and thus the resulting spectroscopic quantities. Applying flexible,^68,69 or even polarizable,^70,71 water models may rectify some of the discrepancy between ODNP experiment and MD simulations.

Additionally, ξ_r is sensitive to the motion of bound waters—such as those buried within soft materials—on longer timescales (100-ps to 10-ns). While we do not expect these longer timescale dynamics to dramatically impact hydration water dynamics in glycerol–water mixtures, the magnitude of the slow time constant τ₁ exceeds 100-ps for x_glyc = 0.3. As the timescales for diffusive and viscous relaxation approach the nanosecond regime, accurately quantifying spectroscopic quantities becomes more difficult for simulations. For glycerol–water, we find that precise quantification of τ₁ becomes challenging for x_glyc > 0.15 [Fig. S7, ESI†].

In the inset of Fig. 3(a), we highlight the spectral density J(ω_S) at the Larmor precession frequency of the electron spin (ω_S = γ_SB₀ ≈ 9.8 GHz) with increasing glycerol concentration. Notably, J(ω_S) is non-monotonic, increasing up to a glycerol concentration of x_glyc = 0.1, then decreasing at higher concentrations [Fig. 3(a) inset]. Accordingly, the computed k_σ values in Fig. 3(c) exhibit the same non-monotonic trend as J(ω_S) with increasing glycerol concentration. Notably, both the experimental and computed values of the relative cross-relaxivity k_σ,r = k_σ/k_σ,pure exhibit the non-monotonic trends with increasing viscosity, with k_σ,r values initially increasing up to a critical glycerol concentration (x_glyc ≈ 0.10) and then decreasing thereafter [Fig. 3(c)]. In the glycerol–water system with compositions from x_glyc = 0 to 0.1, the experimental and MD-derived k_σ,r are in near quantitative agreement. Here, simulation underestimates k_σ,r, but the 95% CI of the MD-derived k_σ,r values bracket the experimental results. For glycerol concentrations beyond x_glyc = 0.1, simulations again underestimate the k_σ,r, values, with the experimental k_σ,r consistently lying above the 95%-CI of the simulation-derived values. We further note a dramatic increase in the relative error in computing k_σ,r at high concentrations.

We attribute the underestimation of k_σ,r to the inability of classical MD simulations to reproduce spectral density amplitudes at high frequency (ω ≫ ω_I). We also note that the agreement between ODNP and MD results is much improved for ξ_r. We believe that this improved agreement results from the dominant influence of spectral density amplitudes at the lower frequencies [e.g., J(ω_I) ≫ J(ω_S)] in the denominator of ξ [eqn (5)]. The agreement of MD-derived values of both ξ_r and k_σ,r with experiments for concentrations between x_glyc = 0 and 0.1 demonstrates that atomistic MD simulations accurately model trends in the translational dynamics of hydration waters on timescales of tens to hundreds of picoseconds.

We also verify the ability of MD simulations to reproduce spectroscopic measures of glycerol–water bulk dynamics via computation of T₁₀[0]. The relative proton longitudinal relaxation time, (T₁₀[0])_r = T₁₀[0]/T₁₀[0]_pure, exhibits a similar monotonic decrease as ξ_r with increasing glycerol content [Fig. 3(d)]. Further, we observe a striking, near-quantitative agreement between experimental and MD-derived (T₁₀[0])_r for x_glyc < 0.10. The MD-derived (T₁₀[0])_r data systematically underestimate the experimental value for x_glyc > 0.075 with the experimental data falling outside of the 95% CI of the MD-calculated values. For x_glyc ≥ 0.1, the computed values underestimate the experimental results by between 20 and 50%. In combination with the ξ_r and k_σ,r results, the (T₁₀[0])_r results suggest that classical atomistic MD simulations can reliably capture the low-frequency contributions to spectroscopic quantities, but less well the high-frequency contribution from instantaneous motion reflected in quantities such as k_σ.

Connecting other hydration water dynamics probes to ODNP measurements

By effectively modeling spectroscopic quantities, we can directly connect ODNP measurements to the microscopic dynamics and structural properties of water. MD simulations also enable computation of other modes of hydration water dynamics that cannot be measured experimentally. For instance, the coupling of these various modes of water dynamics have been studied extensively in liquid water.^32,72–74 One such example is the extended jump model of Laage and Hynes,^27,28 describing the process of water reorientation as dependent simultaneously on the breaking/forming of hydrogen bonds precipitated by large rotational jumps while being rate-limited by translational motion. In this work, we measure characteristic time scales for hydration water translation, rotation, and hydrogen bonding in the glycerol–water mixtures, as summarized in Table S3 (ESI†), using computed autocorrelation functions (ACFs): C⁽⁰⁾_ODNP(t), C_survival(t), C⁽²⁾_OACF(t) and C_HB(t). Further, we explicitly characterize hydration water dynamics via the different relaxation time constants for ODNP-derived diffusion (τ_ODNP), translational diffusion underlying survival probability (τ_survival), rotational diffusion (τ_OACF), and hydrogen bonding (τ_HB).

To complement ODNP-derived measurements of water dynamics, we probe the translational mobility of hydration waters local to the spin probe via the survival probability ACF, C_survival(t) [see Theory Section]. Fig. 4(a) depicts slower decay of C_survival(t) as glycerol concentration increases, indicating the retardation of translational dynamics near the spin probe. We quantify this translational retardation via a characteristic time constant of translational diffusion τ_survival by fitting C_survival(t) to a bi-exponential model c₁e^−t/α₁ + (1 − c₁)e^−t/α₂ and integrating over time. Here, c₁, α₁, and α₂ are the fitting parameters with α₁ > α₁. In agreement with the visible shift in the decay rate of C_survival(t), τ_survival monotonically increases with glycerol concentration [Fig. 4(e)]. There, we illustrate the simultaneous increase in τ_survival and decrease in ξ_r with glycerol concentration. This strong correlation [R² = 0.96] suggests a correspondence between the ODNP coupling factor and measures of the instantaneous translation dynamics of hydration water.

Fig. 4 Molecular dynamics probes of hydration water dynamics correlate strongly with ODNP coupling factors. (a) The survival probability C_survival(t) is the fraction of hydration shell waters that remain continuously within the second hydration shell of the spin probe radical oxygen. (b) The orientational autocorrelation function (OACF) C_OACF(t) measures the rotation of hydration water dipole vectors away from their initial position. (c) The hydrogen bonding survival probability C_HB(t) gives a time scale for water–water hydrogen bond breaking with a hydrogen bond being defined by cutoff radius r_cutoff and cutoff angle θ_cutoff. (d) The ODNP correlation function C_ODNP(t) is used to estimate ODNP spectroscopic quantities. We derive characteristic time constants for (e) translational diffusion τ_survival, (f) rotational diffusion τ_OACF, (g) hydrogen bond lifetimes, and (h) ODNP diffusion τ_ODNP by integrating bi-exponential model fits to the ACFs [(a)–(d), respectively]. Further, these time constants all correlate strongly with relative coupling factor ξ_r.

In addition to the translational mobility, we also examine the hydration water orientational dynamics via the orientational autocorrelation functions (OACFs) [see Theory section]. While we anticipate strong coupling between the translational and rotational mobility of water molecules in dilute solutions, the degree of translation–rotation coupling in crowded glycerol–water mixtures is less intuitive given that the breakdown of orientational–translational coupling has been previously observed for water under confinement.⁷³ As seen in Fig. 4(b), C⁽²⁾_OACF decays more slowly with increasing glycerol content, signifying the retardation of rotational motion of hydration waters. This systematic slowdown of rotational dynamics is further illustrated by the monotonic increase in the characteristic time constants for rotational diffusion τ_OACF. As with τ_survivial we derive τ_OACF by fitting C⁽²⁾_OACF to a bi-exponential model. The simultaneous increase in characteristic timescales τ_survivial and τ_OACF [Fig. 4(e) and (f)] across the whole range of glycerol concentrations establishes that, even in highly crowded glycerol–water mixtures, translational and rotational water diffusion remain strongly coupled. The persistent connection between rotational and translational dynamics suggests that—in contrast to water under geometric confinement (i.e., micelles)—the effect of glycerol on water stems from a distinct physical mechanism such as the development of a collective glycerol–water hydrogen bond network.

Glycerol is often assumed to simply decrease water translational dynamics due to increased mixture viscosity, but the effect of glycerol on water structural dynamics is typically not considered nor well understood. To investigate the relevance of such considerations, we probe the dynamics of hydrogen bonding of spin probe hydration waters via the hydrogen bond survival probability C_HB(t). Like C_survival(t) and C⁽²⁾_OACF(t), we report slower decay of C_HB(t) with rising glycerol content [Fig. 4(c)]. As with τ_survival and τ_OACF, we compute water–water hydrogen bond lifetimes τ_HB by integrating a bi-exponential fit to C_HB(t). The resultant τ_HB values monotonically increase with increasing glycerol concentration, again in a manner strongly correlated with the ODNP results. The simultaneous increase in τ_OACF and τ_HB suggests that glycerol enhances the lifetime of orientationally-coordinated and hydrogen bonded microstructures. Such an enhancement water structure with increasing glycerol content was demonstrated in our recently published study.⁷⁵ In this work, we found that increased glycerol concentration not only slows the water diffusivity, but also enhances water orientational structure by increasing the tetrahedrality of water.⁷⁶ We expand this analysis in the following sub-section, discussing connections between water tetrahedrality and various characteristic time constants in glycerol–water.

With these analyses, we elucidate the as-of-yet unknown persistent coupling between metrics for water translational, rotational, and hydrogen bonding dynamics in glycerol–water, including ODNP-derived and computational translational water diffusivity. We support this finding in Fig. S9 (ESI†) by depicting the connections between the hydration water dynamics relaxation times τ_survival, τ_HB, τ_OACF, and τ_ODNP with spectroscopic quantities k_σ, ξ, and T₁₀[0]. Apart from k_σ, we observe strong correlation between each of these quantities and the others (R² > 0.9). The persistent connection between the characteristic time constants of hydration water dynamics suggests that glycerol, even in significant quantities, does not decouple the translational and rotational modes of water motion, in contrast to dynamic decoupling reported for water under nanoscale confinement⁷³ and in supercooled conditions.^72,74,77 In the context of nanoconfinement between silica planar surfaces, Romero-Vargas Castrillón et al. find that water rotational diffusion is bulk-like near the center of the nanochannel while translational diffusion is suppressed relative to bulk.⁷³ The authors attribute this rotation–translation decoupling to the strong inverse dependences of translational diffusion on density (above bulk near the center channel) and the water–water hydrogen bond network (same as the bulk near the center channel). In contrast, increasing glycerol concentration yields simultaneous increases in mixture density and water–water network hydrogen bond lifetimes (increasing τ_HB), as well as water structuring as measured by tetrahedrality, relative to pure water.

Due to the non-monotonic behavior of k_σ [Fig. 3(c)] with composition, it correlates less well (R² < 0.9) with the monotonically varying relaxation times and spectroscopic quantities. Traditionally, glycerol is thought to act as a classical viscogen without specifically considering its effect on water's molecular structure. With our simulation studies, we observe an increase in the rotational time constant τ_OACF and hydrogen bond lifetimes τ_HB that suggest tetrahedral enhancement with increasing glycerol concentration. To directly quantify tetrahedrality, we require a more direct structural metric, which we pursue next.

Glycerol enhances the population of tetrahedral waters

A number of order parameters have been used to assess the tetrahedrality of water in simulation studies.^78–80 Many such parameters are influenced by the fact that water is undercoordinated near surfaces and in confined environments due to geometrical constraints.^23,78,79,81 For example, the tetrahedral order parameter, q, from Errington and Debenedetti⁸⁰ effectively characterizes water's tetrahedrality in many contexts, but in particular in dilute aqueous mixtures. However, the interpretation of q is more challenging when the average spatial separation between water molecules increases (e.g., concentrated glycerol–water mixtures) because it relies on four nearest water neighbors, which in such cases can populate non-first shell distances such that they adopt less correlated orientational order (decreased tetrahedrality).

As an alternative to q, three-body angle distributions avoid this geometric bias by consider only nearest neighboring water molecules within a distance cutoff of each central water molecule and computing the distribution of neighbor–central–neighbor triplet angles^78,79 averaged over all simulation time steps. We can characterize the population of tetrahedrally coordinated waters, p_tet, by integrating P(θ) over the tetrahedral region of the distribution . In addition to readily accommodating cases when water is under-coordinated, three-body angle distributions have been shown to aptly capture changes in water structure in response to shifts in thermodynamic properties^78,79 and solute chemistry.⁸¹ Furthermore, shifts in three-body angle distribution have also proved to be predictive of solvation thermodynamics for a wide range of colloidal particle sizes.⁷⁸

In Fig. 5(a), we illustrate shifts in the populations of waters in tetrahedral and icosahedral (simple fluid) environments at 109.5 degrees and 64 degrees, respectively, that report on water's structural orientational environment. In Fig. 5(a), we show that the increase in glycerol concentration from x_glyc = 0.01 to 0.3 is accompanied by an increase of 2% in the overall p_tet and an equivalent decrease in the population of icosahedrally-coordinated waters. Curiously, prior work by Monroe and Shell demonstrated that increasing pure water density tended to decrease tetrahedrality as measured by the three body angle distribution,⁷⁸ while glycerol–water behaves in the opposite manner: as mixture density increases (increasing x_glyc), p_tet increases. This finding further supports that glycerol, rather than acting as a simple viscogen, enhances water's structural environment through enhanced water–glycerol hydrogen bonding. Work by Sharp, Vanderkooi, and coworkers supports this finding, using infrared spectroscopy and MD simulation to show that glycerol mimics enhanced hydrogen bonding structures commonly seen in ice.²⁵

Fig. 5 The three-body-angle distribution shows enhanced water tetrahedrality with increasing glycerol concentration. (a) Increasing glycerol concentration in the mixture increases the incidence of tetrahedrally-coordinated waters relative to pure water [P(109.5°) − P_pure(109.5°)] while decreasing the incidence of icosahedrally-coordinated (simple-fluid like) waters [P(64°) − P_pure(64°)]. The increasing population of tetrahedral waters with glycerol concentration correlates strongly to the relative diffusivity of pure water D_H2O/D_H2O,pure at a given mixture composition. Characteristic time constants for (b) translational diffusion, τ_survival, (c) rotational diffusion, τ_OACF, (d) hydrogen bond lifetimes, τ_HB, and (e) the ODNP correlation function, τ_ODNP correlate strongly with R² > 0.99 to the population of tetrahedral waters .

Remarkably, the enhancement of water tetrahedrality strongly correlates with the retardation of hydration water dynamics probes. In Fig. 5(b)–(e), we discover nearly linear relationships (R² > 0.96) between p_tet and several probes of water dynamics, including logarithms of characteristic time constants for translational diffusion τ_survival, rotational diffusion τ_OACF, hydrogen bonding lifetimes τ_HB, and the ODNP correlation time τ_ODNP. The near mono-exponential relationship between p_tet and these time constants demonstrates that structural enhancements—driven by the addition of glycerol—impose systematic retardation of water dynamics. Moreover, the remarkably strong correlation between τ_ODNP and the structural metric p_tet suggests that ODNP measurements indirectly reports on the effect of mixture properties on tetrahedrality of water in bulk solution. The strong correspondence between equilibrium water dynamics and water's molecular structure in glycerol–water is consistent with a previous work by Shell and coworkers.⁷⁶ In this prior work, we find correlation between the system-average water self-diffusivity and several structural metrics such as p_tet, the water–water coordination numbers, water specific volumes, and q among others. Further, we discovered that the system-average water self-diffusivity is accurately predicted by only two structural metrics—including p_tet. In the present work, we expand upon this observation finding strong relationships between not only the system-average translational dynamics of water, but local dynamics for water rotation and water–water hydrogen bonding.

Impact of glycerol on solvation thermodynamics

We apply this understanding of water dynamics and structure to contextualize glycerol–water thermodynamics by characterizing the tendency for a model small hydrophobic molecule (methane) to transfer from an ideal gas to a glycerol–water solution phase in the infinite dilution limit. We directly quantify this by computing the excess free energy of solvation ΔG^ex_solv as discussed in the Methods section. ΔG^ex_solv < 0 suggests favorable solvation of a solute relative to ideal gas phase while the opposite is true for ΔG^ex_solv > 0. In Fig. S10(b) (ESI†), we observe a monotonically increasing and positive ΔG^ex_solv with increasing x_glyc for x_glyc < 0.1. This increase in ΔG^ex_solv is in part explained by the enhancement of water structure [Fig. 5] that presumably generates an increase in the entropic penalty for restructuring water–glycerol. However, for x_glyc > 0.1, ΔG^ex_solv begins to plateau to a constant ΔG^ex_solv ≈ 4.25k_BT. Notably, none of the structural or dynamic metrics above exhibit such a plateau.

To better understand the trend depicted in Fig. S10(b) (ESI†), we apply a free energy decomposition ΔG^ex_solv = U_sw + S_res previously described by Monroe and Shell.⁷⁸ Here, U_sw is the mean interaction energy between methane and glycerol–water in the solvated state (2) and the restructuring entropy S_res is strictly positive and largely gives the entropic penalty to create a void large enough to accommodate the solvation of methane from the gas phase (1) into the solution phase (2). We quantify U_sw by directly calculating the interaction energy between methane and glycerol–water U_sw while S_res follows from S_res = ΔG^ex_solv − U_sw. We observe a monotonic decrease in U_sw with increased glycerol concentration [Fig. 6(a)], indicating a more favorable enthalpy of solvation for methane in glycerol–water. On the other hand, we demonstrate that S_res increases monotonically with glycerol concentration [Fig. 6(b)], which indicates an increase in restructuring penalty to methane solvation. The entropic penalty (S_res > 0) and enthalpic gain (U_sw < 0) for methane solvation in glycerol–water are consistent with the long-known hydration behavior of small hydrophobic molecules in aqueous environments.^82,83 As ΔG^ex_solv becomes less favorable (more positive) for x_glyc < 0.1, U_sw decreases such that further increases in S_res yield constant ΔG^ex_solv for x_glyc > 0.1.

Fig. 6 Decomposition of the solvation free energy of methane into glycerol–water mixtures. We calculate the solvation free energy for a methane molecule via expanded ensemble calculations, decomposing the resultant solvation free energy [Fig. S10(b), ESI†] into (a) enthalpic contribution via the direct energy term U_sw and (b) entropy of solution restructuring S_res. (a) U_sw decreases as more glycerol is added to the mixture. (b) S_res increases as more glycerol is added to the mixture. (c) We observe a strong correlation (R² = 0.99) between the MD-computed ODNP coupling factor ξ_r and S_res.

While the excess solvation free energy, ΔG^ex_solv, does not show a simple correlation with dynamic or structural metrics, the entropic and energetic contributions display striking connections to dynamic metrics. Specifically, we note that the increasing penalty for restructuring glycerol–water S_res displays a strong negative correlation with the MD-computed relative coupling factor ξ_r [R² = 0.99; Fig. 6(c)]. Given that increases in S_res are driven by an enhancement in the underlying glycerol–water solution structure, the connection between S_res and ξ_r again suggests a persistent structure–dynamics connection in glycerol–water. Beyond reiterating the structure–dynamics connections depicted in Fig. 5, the S_res − ξ_r relationship is reminiscent of entropic theories for microscopic dynamics like the viscosity–entropy relationship proposed by Adam and Gibbs.⁸⁴

Adam–Gibbs theory proposes that retardation of liquid dynamics (such as liquid viscosity, η, or self-diffusivity, D) stems from a decrease in the number of available configurational states of a liquid via the configurational entropy S_conf. Notably, a recent study by Handle and Sciorno⁸⁵ demonstrated that a persistent ln D − (TS_conf)⁻¹ connection—per Adam–Gibbs theory—in simulations of pure TIP4P/2005 water. The observed correlation between the probe of water dynamics (ξ_r) and restructuring entropy hints at a theoretical analog to Adam–Gibbs relating equilibrium water dynamics to solvation thermodynamics. If such a S_res-dynamics relationship persists for contexts beyond glycerol–water, it may be possible to forecast the solvation thermodynamics of small molecules in aqueous mixtures without computationally laborious free energy calculations. Probing dynamics–structure–thermodynamics correspondence at heterogeneous surfaces would be fascinating extensions of the present analysis. For instance, one could specifically interrogate how water–glycerol's molecular scale structure and dynamics mediate hydration free energies at protein–water interfaces to elucidate the cosolvent-mediated mechanisms of protein cryoprotection.^86–89

Conclusions

In the present study, we reproduce critical ODNP spectroscopic measures of translational hydration water dynamics using classical atomistic molecular simulations to model equilibrium dynamics in glycerol–water mixtures with semi-quantitative agreement. Further, using the molecular scale detail revealed through MD simulations, we discover strong correlations between ODNP-measured coupling factors and computational probes of translational, rotational, and hydrogen bonding dynamics. The strong relationships found between MD-derived measures of water structure and dynamics are exciting for the potential of ODNP to serve as a surrogate probe of underlying solution structure. Finally, the clear connection between water tetrahedrality, water dynamics, and solvent restructuring entropy suggests a novel framework to describe and quantify the molecular scale mechanisms underlying hydrophobic hydration.

Though the existing literature on local water properties and water-mediated behavior is extensive, our analyses are unique because we directly link simulation to experiment and develop a dynamics–structure–thermodynamics connection. The use of ODNP as a tool to inspect ambient local water environments will greatly complement alternative NMR approaches to studying confined and interfacial water.^90–92 The generality of the experimental and computational methods discussed here will enable further investigations of other systems of broad interest such as water–alcohol and multicomponent mixtures used for protein stabilization (e.g., water–glycerol–DMSO). Further, the dynamics–structure–thermodynamics framework will aide further understanding of the molecular scale mechanisms (e.g., water tetrahedrality and diffusivity) underlying hydration properties in a wide range of chemically and topologically heterogeneous interfaces such as at protein–water interfaces.

Author contributions

Dennis Robinson Brown (conceptualization, MD simulations and methodology, data analysis and analytical modelling, writing – original draft); Thomas R. Webber (conceptualization, ODNP experiments and methodology, data analysis, writing – original draft); Thomas M. Casey (conceptualization, ODNP experiments and methodology, data analysis, writing – review & editing); John Franck (conceptualization, analytical modelling, writing – review & editing); M. Scott Shell (conceptualization, writing – review & editing, supervision, funding acquisition); Songi Han (conceptualization, writing – review & editing, supervision, funding acquisition).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported as part of the Center for Materials for Water and Energy Systems (M-WET), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award #DE-SC0019272. Use was made of computational facilities purchased with funds from the National Science Foundation (CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF DMR 1720256) at UC Santa Barbara. Support for ODNP method development comes from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC2033 390677874 RESOLV. Support for surface hydration structure and dynamics is supported by the NIH MIRA grant R35GM136411 awarded to SH.

References

1. N. B. Rego, E. Xi and A. J. Patel, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2018234118.
2. S. N. Jamadagni, R. Godawat and S. Garde, Langmuir, 2009, 25, 13092–13099.
3. R. Godawat, S. N. Jamadagni and S. Garde, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15119–15124.
4. J. I. Monroe, S. Jiao, R. J. Davis, D. Robinson Brown, L. E. Katz and M. S. Shell, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2020205118.
5. T. Young, R. Abel, B. Kim, B. J. Berne and R. A. Friesner, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 808–813.
6. H. Wirtz, S. Schäfer, C. Hoberg and M. Havenith, J. Infrared Milli Terahz Waves, 2018, 39, 816–827.
7. I. Moskowitz, M. A. Snyder and J. Mittal, J. Chem. Phys., 2014, 141, 18C532.
8. D. Laage, G. Stirnemann and J. T. Hynes, J. Phys. Chem. B, 2009, 113, 2428–2435.
9. J. I. Monroe and M. S. Shell, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 8093–8098.
10. K. Lum, D. Chandler and J. D. Weeks, J. Phys. Chem. B, 1999, 103, 4570–4577.
11. J. Mittal and G. Hummer, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 20130–20135.
12. S. Garde, G. Hummer, A. E. García, L. R. Pratt and M. E. Paulaitis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1996, 53, R4310–R4313.
13. H. Acharya, S. Vembanur, S. N. Jamadagni and S. Garde, Faraday Discuss., 2010, 146, 353–365.
14. E. Xi, V. Venkateshwaran, L. Li, N. Rego, A. J. Patel and S. Garde, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 13345–13350.
15. N. Giovambattista, P. G. Debenedetti and P. J. Rossky, J. Phys. Chem. C, 2007, 111, 1323–1332.
16. R. Barnes, S. Sun, Y. Fichou, F. W. Dahlquist, M. Heyden and S. Han, J. Am. Chem. Soc., 2017, 139, 17890–17901.
17. G. Schirò, Y. Fichou, F.-X. Gallat, K. Wood, F. Gabel, M. Moulin, M. Härtlein, M. Heyden, J.-P. Colletier, A. Orecchini, A. Paciaroni, J. Wuttke, D. J. Tobias and M. Weik, Nat. Commun., 2015, 6, 6490.
18. J. G. Davis, K. P. Gierszal, P. Wang and D. Ben-Amotz, Nature, 2012, 491, 582–585.
19. F. Böhm, G. Schwaab and M. Havenith, Angew. Chem., Int. Ed., 2017, 56, 9981–9985.
20. N. B. Rego, E. Xi and A. J. Patel, J. Am. Chem. Soc., 2019, 141, 2080–2086.
21. A. M. Schrader, J. I. Monroe, R. Sheil, H. A. Dobbs, T. J. Keller, Y. Li, S. Jain, M. S. Shell, J. N. Israelachvili and S. Han, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 2890–2895.
22. J. M. Franck, Y. Ding, K. Stone, P. Z. Qin and S. Han, J. Am. Chem. Soc., 2015, 137, 12013–12023.
23. J. Monroe, M. Barry, A. DeStefano, P. Aydogan Gokturk, S. Jiao, D. Robinson-Brown, T. Webber, E. J. Crumlin, S. Han and M. S. Shell, Ann. Rev. Chem. Biomol. Eng., 2020, 11, 523–557.
24. J. M. Franck and S. Han, in Methods in Enzymology, ed. A. J. Wand, Academic Press, 2019, vol. 615, pp. 131–175.
25. J. L. Dashnau, N. V. Nucci, K. A. Sharp and J. M. Vanderkooi, J. Phys. Chem. B, 2006, 110, 13670–13677.
26. M. Heyden, J. Chem. Phys., 2019, 150, 094701.
27. D. Laage, Science, 2006, 311, 832–835.
28. D. Laage and J. T. Hynes, J. Phys. Chem. B, 2008, 112, 14230–14242.
29. S. Xiao, F. Figge, G. Stirnemann, D. Laage and J. A. McGuire, J. Am. Chem. Soc., 2016, 138, 5551–5560.
30. H. F. M. C. Martiniano and N. Galamba, J. Phys. Chem. B, 2013, 117, 16188–16195.
31. G. Stirnemann, S. R.-V. Castrillón, J. T. Hynes, P. J. Rossky, P. G. Debenedetti and D. Laage, Phys. Chem. Chem. Phys., 2011, 13, 19911.
32. P. Tan, J. Huang, E. Mamontov, V. García Sakai, F. Merzel, Z. Liu, Y. Ye and L. Hong, Phys. Chem. Chem. Phys., 2020, 22, 18132–18140.
33. B. Qiao, F. Jiménez-Ángeles, T. D. Nguyen and M. Olvera de la Cruz, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 19274–19281.
34. J. N. Dahanayake and K. R. Mitchell-Koch, Phys. Chem. Chem. Phys., 2018, 20, 14765–14777.
35. S.-H. Chen, P. Gallo, F. Sciortino and P. Tartaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1997, 56, 4231–4243.
36. J.-C. Perrin, S. Lyonnard and F. Volino, J. Phys. Chem. C, 2007, 111, 3393–3404.
37. B. D. Armstrong and S. Han, J. Am. Chem. Soc., 2009, 131, 4641–4647.
38. Y. Xu and M. Havenith, J. Chem. Phys., 2015, 143, 170901.
39. J. Zhang, J. Phys. Chem. B, 2019, 123, 2971–2977.
40. H. Hoshina, Y. Iwasaki, E. Katahira, M. Okamoto and C. Otani, Polymer, 2018, 148, 49–60.
41. S. J. Kim, B. Born, M. Havenith and M. Gruebele, Angew. Chem., Int. Ed., 2008, 47, 6486–6489.
42. C. P. Lawrence and J. L. Skinner, J. Chem. Phys., 2003, 118, 264–272.
43. H. Moon, R. P. Collanton, J. I. Monroe, T. M. Casey, M. S. Shell, S. Han and S. L. Scott, J. Am. Chem. Soc., 2022, 144, 1766–1777.
44. D. D. Mahanta, D. R. Brown, S. Pezzotti, S. Han, G. Schwaab, M. S. Shell and M. Havenith, Chem. Sci., 2023, 14, 7381.
45. I. Kaminker, R. Barnes and S. Han, in Methods in Enzymology, ed. P. Z. Qin and K. Warncke, Academic Press, 2015, vol. 564, pp. 457–483.
46. J. M. Franck, A. Pavlova, J. A. Scott and S. Han, Prog. Nucl. Magn. Reson. Spectrosc., 2013, 74, 33–56.
47. S. Izadi, R. Anandakrishnan and A. V. Onufriev, J. Phys. Chem. Lett., 2014, 5, 3863–3871.
48. J. Blieck, F. Affouard, P. Bordat, A. Lerbret and M. Descamps, Chem. Phys., 2005, 317, 253–257.
49. R. Chelli, P. Procacci, G. Cardini and S. Califano, Phys. Chem. Chem. Phys., 1999, 1, 879–885.
50. D. A. Jahn, F. O. Akinkunmi and N. Giovambattista, J. Phys. Chem. B, 2014, 118, 11284–11294.
51. D. A. Case, 2018.
52. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, Williams, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, 2016.
53. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
54. Development of the Second Generation of the General AMBER Force Field, https://www.researchgate.net/project/Development-of-the-Second-Generation-of-the-General-AMBER-Force-Field, (accessed June 17, 2020).
55. D. Sezer, J. H. Freed and B. Roux, J. Phys. Chem. B, 2008, 112, 5755–5767.
56. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
57. P. Eastman and V. Pande, Comput. Sci. Eng., 2010, 12, 34–39.
58. M. R. Shirts and J. D. Chodera, J. Chem. Phys., 2008, 129, 124105.
59. D. Sezer, M. Gafurov, M. J. Prandolini, V. P. Denysenkov and T. F. Prisner, Phys. Chem. Chem. Phys., 2009, 11, 6638.
60. L.-P. Hwang and J. H. Freed, J. Chem. Phys., 1975, 63, 4017.
61. J. H. Ortony, B. Qiao, C. J. Newcomb, T. J. Keller, L. C. Palmer, E. Deiss-Yehiely, M. Olvera de la Cruz, S. Han and S. I. Stupp, J. Am. Chem. Soc., 2017, 139, 8915–8921.
62. D. Sezer, M. J. Prandolini and T. F. Prisner, Phys. Chem. Chem. Phys., 2009, 11, 6626.
63. B. D. Armstrong, P. Soto, J.-E. Shea and S. Han, J. Magn. Reson., 2009, 200, 137–141.
64. B. D. Armstrong and S. Han, J. Chem. Phys., 2007, 127, 104508.
65. N. Bloembergen, E. M. Purcell and R. V. Pound, Phys. Rev., 1948, 73, 679–712.
66. S. Shin and A. P. Willard, J. Chem. Theory Comput., 2018, 14, 461–465.
67. A. Luzar and D. Chandler, J. Chem. Phys., 1993, 98, 8160–8173.
68. Y. Wu, H. L. Tepper and G. A. Voth, J. Chem. Phys., 2006, 124, 024503.
69. M. A. González and J. L. F. Abascal, J. Chem. Phys., 2011, 135, 224516.
70. T. E. Gartner, K. M. Hunter, E. Lambros, A. Caruso, M. Riera, G. R. Medders, A. Z. Panagiotopoulos, P. G. Debenedetti and F. Paesani, J. Phys. Chem. Lett., 2022, 13, 3652–3658.
71. V. Babin, C. Leforestier and F. Paesani, J. Chem. Theory Comput., 2013, 9, 5395–5403.
72. S. Dueby, V. Dubey and S. Daschakraborty, J. Phys. Chem. B, 2019, 123(33), 7178–7189.
73. S. Romero-Vargas Castrillón, N. Giovambattista, I. A. Aksay and P. G. Debenedetti, J. Phys. Chem. B, 2009, 113, 1438–1446.
74. V. Dubey, S. Erimban, S. Indra and S. Daschakraborty, J. Phys. Chem. B, 2019, 123, 10089–10099.
75. D. C. Robinson Brown, T. R. Webber, S. Jiao, D. M. Rivera Mirabal, S. Han and M. S. Shell, J. Phys. Chem. B, 2023, 127(20), 4577–4594.
76. D. C. Robinson Brown, T. R. Webber, S. Jiao, D. M. Rivera Mirabal, S. Han and M. S. Shell, J. Phys. Chem. B, 2023, 127, 4577–4594.
77. P. G. Debenedetti and F. H. Stillinger, Nature, 2001, 410, 259–267.
78. J. I. Monroe and M. S. Shell, J. Chem. Phys., 2019, 151, 094501.
79. A. Chaimovich and M. S. Shell, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81(6), 060104.
80. J. R. Errington and P. G. Debenedetti, Nature, 2001, 409, 318–321.
81. P. Stock, J. I. Monroe, T. Utzig, D. J. Smith, M. S. Shell and M. Valtiner, ACS Nano, 2017, 11, 2586–2597.
82. E. Gallicchio, M. M. Kubo and R. M. Levy, J. Phys. Chem. B, 2000, 104, 6271–6285.
83. H. S. Frank and M. W. Evans, J. Chem. Phys., 1945, 13, 507–532.
84. G. Adam and J. H. Gibbs, J. Chem. Phys., 1965, 43, 139–146.
85. P. H. Handle and F. Sciortino, Mol. Phys., 2018, 116, 3366–3371.
86. J. J. Towey, A. K. Soper and L. Dougan, Faraday Discuss., 2013, 167, 159.
87. L. Weng, C. Chen, J. Zuo and W. Li, J. Phys. Chem. A, 2011, 115, 4729–4737.
88. L. Weng, S. L. Stott and M. Toner, Annu. Rev. Biomed. Eng., 2019, 21, 1–31.
89. J. J. Towey, A. K. Soper and L. Dougan, J. Phys. Chem. B, 2012, 116, 13898–13904.
90. F. Persson, P. Söderhjelm and B. Halle, J. Chem. Phys., 2018, 148, 215103.
91. N. V. Nucci, M. S. Pometun and A. J. Wand, Nat. Struct. Mol. Biol., 2011, 18, 245–249.
92. C. Mattea, J. Qvist and B. Halle, Biophys. J., 2008, 95, 2951–2963.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp00030g

---