The properties of aqueous solutions of ionic, zwitterionic, and polar solutes are of interest to many fields. For instance, one of the many anomalous properties of aqueous solutions is the behavior of water diffusion in different monovalent salt solutions. In addition, solutes can affect the stabilities of macromolecules such as proteins in aqueous solution. Here, the diffusivities of aqueous solutions of sodium chloride, potassium chloride, tri-methylamine oxide (TMAO), urea, and TMAO-urea are examined in molecular dynamics simulations. The decrease in the diffusivity of water with the concentration of simple ions and urea can be described by a simple model in which the water molecules hydrogen bonded to the solutes are considered to diffuse at the same rate as the solutes, while the remainder of the water molecules are considered to be bulk and diffuse at almost the same rate as pure water. On the other hand, the decrease in the diffusivity of water with the concentration of TMAO is apparently affected by a decrease in the diffusion rate of the bulk water molecules in addition to the decrease due to the water molecules hydrogen bonded to TMAO. In other words, TMAO enhances the viscosity of water, while urea barely affects it. Overall, this separation of water molecules into those that are hydrogen bonded to solute and those that are bulk can provide a useful means of understanding the short- and long-range effects of solutes on water.

Water and aqueous solutions display numerous anomalous and unusual properties that distinguish water from other liquids. One of the anomalous behaviors of water in aqueous salt solutions is that the concentration dependence of the diffusivity appears to be a function of ion size. Diffusion coefficients have been determined for water in electrolyte solutions by proton magnetic relaxation rates1 and quasielastic neutron scattering spectroscopy2 and for both ions and water in electrolyte solutions by isotopic labeling.3 These experiments show that while the diffusion of water decreases with concentration in small salts up to 4 to 5M salt, it increases first before it decreases in large salts.1,3 Interestingly, the change in behavior appears to occur between NaCl and KCl solutions2 which are also important electrolytes in biology. Understanding this phenomena has been hindered because additive and polarizable empirical models in conventional molecular dynamics (MD) simulations predict a decrease with salt concentration, regardless of the cation,4 although ab initio MD (AIMD) simulations are able to predict this effect.5 More recent work has attributed this effect to the charge transfer between ions and water6,7 using the discrete charge-transfer (DCT) method8,9 and to background electronic polarization10 using a scaled-ionic-charge model.

In addition, a variety of small solutes in aqueous solutions have been shown to exhibit both stabilizing and destabilizing effects on the water structure, the so-called kosmotropic and chaotropic effects, respectively. These solutes also tend to stabilize and destabilize, respectively, macromolecules such as proteins. Since both ions11,12 and polar13 molecules can be either kosmotropic or chaotropic, understanding what leads to their effects is the subject of many studies. For instance, tri-methylamine oxide (TMAO) stabilizes proteins, while urea destabilizes proteins. Thus, simulations of TMAO14,15 and urea16,17 in water have been used to gain a molecular understanding of the effects of each on water structure, although the effects of urea on the water structure is still under debate.18–20 In addition, TMAO counteracts the effects of urea on protein denaturation at a 1:2 ratio of concentrations of TMAO to urea.21,22 Interestingly, sharks, rays, and other teleosts may use 1:2 mixture of TMAO and urea to protect against high hydrostatic pressure, since they have this mixture at concentrations directly proportional to the depth at which they are found.23,24 Simulations have also been used to examine the counteracting effects of TMAO on urea,25,26 although the cause is still not clear.

Here, the focus is on understanding how different types of solutes can influence the dynamics and structure of water. The diffusion coefficients in two ionic aqueous solutions, NaCl and KCl, and binary and ternary mixtures of a zwitterion and a polar solute, namely, TMAO and urea, in water are examined in MD simulations. While comparisons are made to the available experiment for reference, the main purpose here is to present a framework to understand how different types of solutes can affect hydration and bulk water, rather than the accuracy of the force field. The two ionic aqueous solutions are presented first, followed by the TMAO and urea aqueous solutions. For each, the decrease in the diffusivity of water is related to the number of water molecules hydrogen bonded to the solutes, which are considered to diffuse at the same rate as the solute, and the remainder of the water molecules, or bulk water, which are considered to diffuse at a separate rate. These simulations point to differences in the behavior of bulk water for the different solutions.

Molecular dynamics simulations were performed using the molecular mechanics package CHARMM c41a2.27 The procedure was mostly the default protocol in CHARMMing;28 pertinent details including changes from the default protocols are noted here. Water was modeled using TIP4P-Ew.29 Sodium, potassium, and chloride ions were represented using the parameters optimized for TIP4P-Ew.30 TMAO and urea were represented using the CHARMM36 all-atom force field.31 The SHAKE algorithm32 was used to fix the length of the covalent bonds for hydrogen atoms. The ranges for nonbonded interactions were somewhat different for the ion simulations and the TMAO/urea simulations; the latter used the CHARMM defaults and so appear in parentheses here. For nonbonded interactions, pair and image lists were generated out to 12 Å (16 Å) and updated heuristically. The van der Waals interactions were turned off using the CHARMM switching function from 8 to 10 Å (10 to 12 Å). The electrostatic interactions were calculated using Particle Mesh Ewald (PME) with a β-spline coefficient of 6, a κ value of 0.43 Å−1 (0.34 Å−1), and 108 (80) k-vectors in each Cartesian direction.

Cubic simulation boxes were prepared from equilibrated boxes of pure water. For the ionic solutions, the original box had 4096 water molecules and a side length of ∼50 Å at concentrations of about 0.15, 1, 2, 3, and 4M. For the TMAO and urea solutions, the original box had ∼2280 water molecules and a side length of ∼40 Å at concentrations of about 0.5, 1, 1.5, 2, 3, 4M plus additional higher concentrations. The numbers of each type of molecules and exact concentrations in the systems are given in Tables S1–S5 of the supplementary material. In addition, each TMAO/urea system was subjected to initial energy minimization (200 steps of the steepest descent method). During the energy minimization, heavy atoms of water and co-solute molecules were harmonically restrained with a spring constant of 2 kcal/(mol-Å2). This constraint was removed afterwards.

The simulations employed the leapfrog Verlet algorithm33 with a time step of 1 fs, and NVT and NPT ensembles were maintained using the Nosé-Hoover algorithm.34,35 For ions, after heating from 98 K to 298 K in 5 K increments for 100 ps in the NVT ensemble, the systems were allowed to run unperturbed in the NPT ensemble at a constant pressure of 1 atm and a temperature of 298 K for 0.5 ns of equilibration followed by 5.5 ns of production run. For TMAO/urea systems, after heating from 0 K to 300 K in 30 ps, the system was scaled every 0.2 ps for 170 ps, in the NPT ensemble. Next, the systems were allowed to run unperturbed in the NPT ensemble at a constant pressure of 1 atm and a temperature of 300 K for 5 ns of equilibration followed by 25 to 40 ns of production run (see Tables S3–S5 of the supplementary material for details). For all simulations, coordinates were saved every 1 ps for analysis.

The water molecules considered to be hydrogen bonded to a solute used a distance cutoff of 3.45 Å for Na+–O, 3.80 Å for K+–O, 2.4 Å for O–H or N–H, and 3.0 Å for Cl–H, with no angle cutoff. Thus, water molecules in the first hydration shell of cations are referred to as “hydrogen bonded” for brevity, although strictly it is an electrostatic interaction. The average lifetime of a hydrogen bond is defined as the average duration that a pair of molecules continuously satisfy the criteria for being hydrogen bonded. The limits of the integration of the radial distribution to determine the coordination numbers of the ion pairs were 3.85 Å for K+–Cl and 3.50 Å for Na+–Cl.

The average diffusion coefficient of all water in each simulation was calculated. In addition, diffusion coefficients for two subpopulations of water in the simulations were calculated: DhydMD is the average diffusion coefficient of all water molecules that are considered to be hydrogen bonded to the solute and DbulkMD is the average diffusion coefficient of all of the rest of the water molecules. As in previous work,36 diffusion coefficients were corrected for box size using the method of Yeh and Hummer37 except that experimental viscosities for NaCl,38 KCl,39 TMAO,40 and urea41 in aqueous solution were used. For TMAO/urea, the diffusion coefficient was calculated for every 5-ns interval to obtain error bars.

The diffusion coefficients are analyzed in terms of a simple “hydrogen bonded-bulk” (HB-B) model for the diffusion of water in solutions. First, water that is hydrogen bonded to a solute is assumed to diffuse at the same rate as the solute, while the rest of the water is assumed to diffuse as bulk water. Thus, for a solute i,

Dwc=1NwciNicniwcDic+11NwciNicniw(c)Dbulkc,
(1)

where c is the solute concentration, Ni is the number of i, Di is the diffusion coefficient of i, and niw is the number of water molecules hydrogen bonded to i. In addition, the diffusion coefficient of bulk water at a given concentration, Dbulk(c), is given by

Dbulk(c)=Dw(0)1iaici,
(2)

where Dw(0) is the diffusion coefficient for pure water, ci is the concentration of the solute i, and ai is a coefficient for the concentration dependence of i. The ai coefficients are obtained from the percentage decrease in the diffusivity of i with the concentration from the relevant binary mixture of i and water by assuming that values of Di(0) and Dbulk(0) are controlled by the diffusing object, while the concentration dependences are controlled by the viscosity of the media. In other words, the ai coefficients are obtained by fitting the diffusion coefficients of the solute to the following equation:

Di(c)=Di(0)1aici,
(3)

for the binary solution of solute i and water. For the cases studied here, niw(c), the number of hydrogen-bonded waters for the solutes, can be assumed to be the values at infinite dilution niw(0), regardless of solute concentration, since relatively little concentration dependence is observed. In addition, for most of the cases studied here, Dbulk(c), the diffusion coefficient of bulk water, can be assumed to be that for pure water Dw(0) and Di(c), the diffusion coefficient of the solute, can be assumed to be the values at infinite dilution Di(0).

The results are presented first for NaCl and KCl in TIP4P-Ew water and then for binary and ternary solutions of TMAO and urea in TIP4P-Ew water.

The number of hydrogen-bonded water molecules to each ion were calculated at different concentrations and compared with the experimental results (Fig. 1). For brevity, water molecules that are in the first hydration shell of cations are referred to as being hydrogen-bonded. The experimental values had been obtained by integrating radial distribution functions from neutron diffraction data with hydrogen isotope substitution by an analysis using the empirical potential structure refinement technique and so are actually coordination numbers.42 As a whole, the number of hydrogen-bonded waters for ions in TIP4P-Ew appears to be too high compared to the values from experiment. In particular, although the experimental bars are large, the simulation values for Na+ and K+ (5.89 and 6.5, respectively, at 1M) and for Cl in both NaCl and KCl solutions (6.7 and 6.5, respectively, at 1M) are consistently larger by over ∼0.5 water molecule, barely falling within the upper limit of the error bars. However, the coordination numbers of water around Cl in both NaCl and KCl solutions are close to the experimental values (see Fig. S1 in the supplementary material) which are actually coordination numbers as mentioned above. Thus, the structure of the hydration shell of Cl in the simulations appears correct and some water molecules that are considered hydrogen bonded to Cl according to the criteria in Sec. II actually fall in the second coordination shell. However, since the criteria for “hydrogen-bonded” water molecules to ions is actually the same as the criteria for being coordinated, the number of waters around the cations in the simulations appears to be too large. As the concentration increases, the salts tend to cluster, so fewer water molecules can hydrogen bond to the ions. Since the coordination numbers of Cl increase to 0.17 around Na+ and 0.63 around K+ at 4M, the number of hydrogen bonded waters around the ions decreases more with ion concentration in KCl solutions.

FIG. 1.

Number of hydrogen-bonded water molecules to cations (squares) and anions (circles) as a function of ionic concentration: (a) NaCl and (b) KCl aqueous solutions from molecular dynamics simulations (green) and from experiment42 (black).

FIG. 1.

Number of hydrogen-bonded water molecules to cations (squares) and anions (circles) as a function of ionic concentration: (a) NaCl and (b) KCl aqueous solutions from molecular dynamics simulations (green) and from experiment42 (black).

Close modal

Next, the diffusivities in both NaCl and KCl solutions as a function of concentration were compared with experiment. As a whole, the diffusion coefficients for the ions, Dion, in TIP4P-Ew are too slow by ∼0.5 × 10−9 m2/s compared to experiment (Fig. 2). This may be because the simulation overpredicts the number of water molecules hydrogen bonded to the ions, which increases the size of the diffusing body and lowers its diffusion coefficient. The diffusion coefficients for water Dw are generally closer to experiment, but the value for pure TIP4P-Ew water is too high compared to experiment (2.61 vs. 2.3 × 10−9 m2/s, respectively) and the slope of the decrease with concentration is too large (Fig. 3). Moreover, the almost concentration independence of Dw for KCl solutions is not seen in the simulation.

FIG. 2.

Diffusion coefficients as a function of ionic concentration of cations (squares) and anions (circles) in (a) NaCl and (b) KCl aqueous solutions from molecular dynamics simulations (green) and from experiment3 (black). The estimated experimental precision was 0.5% except for the data for Na+, which was a fit to several data sets with a standard error of 0.6%.

FIG. 2.

Diffusion coefficients as a function of ionic concentration of cations (squares) and anions (circles) in (a) NaCl and (b) KCl aqueous solutions from molecular dynamics simulations (green) and from experiment3 (black). The estimated experimental precision was 0.5% except for the data for Na+, which was a fit to several data sets with a standard error of 0.6%.

Close modal
FIG. 3.

Diffusion coefficients of water as a function of ionic concentration in (a) NaCl and (b) KCl aqueous solutions from molecular dynamics simulations (green solid symbol) and from experiment1,2 (black solid symbols). Experimental error was ±2%. Results are also shown using the HB-B model using data from simulations (green open symbols) and experiment (black open symbols).

FIG. 3.

Diffusion coefficients of water as a function of ionic concentration in (a) NaCl and (b) KCl aqueous solutions from molecular dynamics simulations (green solid symbol) and from experiment1,2 (black solid symbols). Experimental error was ±2%. Results are also shown using the HB-B model using data from simulations (green open symbols) and experiment (black open symbols).

Close modal

The diffusion coefficients of water in the salt solutions are examined in terms of the simple HB-B model [Eq. (1)] where the diffusion coefficient of bulk water Dbulk(c) is assumed to be the same as that for pure water Dw(0), regardless of ionic concentration. Thus, for the solutes i = Na+, K+, and Cl, ai are assumed to be zero in Eqs. (2) and (3). The Dw(c) predicted by the simple HB-B model using Di and niw from either simulations or experiments (values at infinite dilution given in Table S6 of the supplementary material) correlate remarkably well with the respective simulation or experimental Dw(c) (Fig. 3), which indicates that the HB-B model is a reasonable framework to understand the average diffusivity of water. This is especially remarkable for the experiment, given the large error bars for the coordination numbers. The HB-B model somewhat underestimates the decrease in Dw(c) for NaCl solutions seen in simulations and does not capture the initial increase in Dw(c) for KCl solutions seen in experiment, which will be discussed further in Sec. IV.

The average number of hydrogen bonds between the solutes, TMAO and urea, and water were calculated for aqueous solutions of TMAO, urea, and a 1:2 ratio of TMAO-urea; to our knowledge, no experimental numbers exist for either hydrogen bond numbers or coordination numbers. The oxygen of TMAO can accept 3.4 hydrogen bonds from water at 0.5M. The average number of hydrogen-bonded water molecules to TMAO is only slightly dependent on concentration in the TMAO solution, while the number decreases slightly in the TMAO-urea solutions, which can be attributed to TMAO hydrogen bonding to urea instead of water as indicated by the increasing number of TMAO-urea hydrogen bonds [Fig. 4(a)]. Urea can accept 2.3 hydrogen bonds to water at 0.5M via its oxygen or donate a total of 2.5 hydrogen bonds to water at 0.5M via its four nitrogens or it can hydrogen bond to itself. The average number of hydrogen-bonded water molecules to urea has a greater dependence on the total solute concentration mainly because the probability of urea-urea hydrogen bonding increases with urea concentration [Fig. 4(b)]. The effects of urea hydrogen bonding to TMAO instead of water are much less (<∼0.2 at any concentration). In addition, when the results are plotted as a function of the total solute concentration, the number of water molecules hydrogen bonded to TMAO is almost the same in either the TMAO or the TMAO-urea aqueous solution and the number water molecules hydrogen bonded to urea is almost the same in either the urea or the TMAO-urea solution, since the same number of water molecules are available for hydrogen bonding at a given total solute concentration (see Fig. S2 in the supplementary material).

FIG. 4.

Number of hydrogen bonds in (a) TMAO (red) and TMAO-urea (magenta) aqueous solutions versus TMAO concentration and (b) urea (blue) and TMAO-urea (magenta) aqueous solutions versus urea concentration from simulations. In key, two letters (where T is TMAO, U is urea, W is water, and for urea, the subscript O indicates that urea is the donor and the subscript N indicates that urea is the acceptor) indicate a hydrogen bond pair with solution indicated in parentheses. Hydrogen bonds to water (squares), urea (circles), and TMAO (triangles) as an acceptor (solid symbols) and as a donor (open symbols) are indicated.

FIG. 4.

Number of hydrogen bonds in (a) TMAO (red) and TMAO-urea (magenta) aqueous solutions versus TMAO concentration and (b) urea (blue) and TMAO-urea (magenta) aqueous solutions versus urea concentration from simulations. In key, two letters (where T is TMAO, U is urea, W is water, and for urea, the subscript O indicates that urea is the donor and the subscript N indicates that urea is the acceptor) indicate a hydrogen bond pair with solution indicated in parentheses. Hydrogen bonds to water (squares), urea (circles), and TMAO (triangles) as an acceptor (solid symbols) and as a donor (open symbols) are indicated.

Close modal

Next, the diffusivities of TMAO and urea in all of the solutions are examined as a function of the concentration of either TMAO or urea, respectively (Fig. 5). The MD results for TMAO in aqueous solutions agree very well with experiment,43 while the MD results for urea are somewhat too fast compared with experiment44 but have a similar decrease with concentration. Perhaps the most remarkable difference of the behavior compared to salts is the large decrease in the diffusion coefficient of TMAO with concentration over a similar range of concentrations. In addition, while replacing water with urea seems to have a rather small effect on further slowing the diffusion of TMAO [Fig. 5(a)], replacing water with TMAO greatly slows down the diffusion of urea [Fig. 5(b)]. The results are also plotted versus the total solute concentration in Fig. S3 of the supplementary material.

FIG. 5.

Diffusion coefficients of (a) TMAO in TMAO (red) and TMAO-urea (magenta) aqueous solutions as a function of TMAO concentration and (b) urea in urea (blue) and TMAO-urea (magenta) aqueous solutions as a function of urea concentration from molecular dynamics simulations. Experimental results are also shown in black for TMAO43 and urea.44 

FIG. 5.

Diffusion coefficients of (a) TMAO in TMAO (red) and TMAO-urea (magenta) aqueous solutions as a function of TMAO concentration and (b) urea in urea (blue) and TMAO-urea (magenta) aqueous solutions as a function of urea concentration from molecular dynamics simulations. Experimental results are also shown in black for TMAO43 and urea.44 

Close modal

Finally, the diffusivities of water in all of the solutions were calculated as a function of total solute concentration (Fig. 6). As before, the diffusion coefficient for pure water calculated for TIP4P-Ew is higher than experimental results; however, the much sharper decrease with concentration for TMAO than urea seen in experiment is predicted well by TIP4P-Ew. While no experimental data exist to our knowledge for the 1:2 mixture of TMAO-urea in aqueous solution, the simulation predicts that it is intermediary between the TMAO in aqueous solution and urea in aqueous solution.

FIG. 6.

Diffusion coefficients of water as a function of total solute concentration in aqueous solutions of TMAO (red), urea (blue), and TMAO-urea (magenta) from molecular dynamics simulations (solid symbols) and from the HB-B model (open symbols and dashed lines). Experimental results are also shown for TMAO (black solid symbols)43 and urea (black open symbols).44 Results are also shown using the HB-B model (open symbols).

FIG. 6.

Diffusion coefficients of water as a function of total solute concentration in aqueous solutions of TMAO (red), urea (blue), and TMAO-urea (magenta) from molecular dynamics simulations (solid symbols) and from the HB-B model (open symbols and dashed lines). Experimental results are also shown for TMAO (black solid symbols)43 and urea (black open symbols).44 Results are also shown using the HB-B model (open symbols).

Close modal

The diffusion coefficients of water in the TMAO and urea solutions are also examined in terms of the HB-B model [Eq. (1)] where the solutes i = TMAO and urea. While the assumption that Dbulk(c) is equal to Dw(0) regardless of solute concentration works fairly well for urea, it is apparent that a decrease in Dbulk(c) with solute concentration must be invoked for TMAO in Eq. (2) so that aTMAO is nonzero. Thus, Di, niw, and aTMAO were obtained from either simulations or experiments as appropriate (values at infinite dilution given in Table S7 of the supplementary material). In addition, the parameters for a solute i were obtained from the binary solutions of solute i and water, and the same values were used for the tertiary solutions. The Dw(c) predicted by the HB-B model correlate remarkably well with the respective simulation or experimental Dw(c) (Fig. 6). Since the estimated amount of bulk water is approximately one half of all water at ∼3M total solute concentration, the agreement is remarkable even at high concentrations.

Overall, the agreement between the Dw(c) predicted by the HB-B model and that calculated directly from the simulations of ions, TMAO, and urea in aqueous solution suggests that Dw(c) can be considered as contributions from water hydrogen bonded to the solutes that diffuse at the same rate as the solute and from bulk water. Of course, the simulation results are dependent on the force field, and given known discrepancies especially for force field descriptions of ions in water, the interpretations of simulation results presented here are not meant as definitive interpretations of what happens in the real system. However, this analysis can lead to understanding where disagreement between simulations and experiment may lie, which can aid force field development.

One assumption for the HB-B model is the number of water molecules diffusing with the solute, which is assumed to be the number that are hydrogen bonded to it. This can be examined in light of hydrogen bond lifetimes calculated from the simulation (see Sec. II). For the salt solutions, since the ion-water hydrogen bond lifetimes are ∼37 ps for Na+ and ∼7 ps for K+ and Cl compared to the water-water hydrogen bond lifetime of 3 ps for TIP4P-Ew, the assumption that water hydrogen bonded to an ion diffuses at the same rate as the ion appears reasonable. For TMAO solutions, since the average TMAO-water hydrogen bond lifetime is 12–30 ps as concentration increases to almost 6M TMAO compared to water-water hydrogen bond lifetime of 3 ps for TIP4P-Ew, the assumption that hydrogen bonded water diffuses at the same rate as the solute appears reasonable. However, for urea solutions, the average urea-water hydrogen bond lifetime is 3.6–4.6 ps as an acceptor and ∼1.5 ps as a donor, so the water molecules hydrogen bonded to especially the NH might not stay long enough to affect the net diffusion rate of urea as much. In addition, the assumption of the concentration independence of the number of water molecules diffusing with the solute appears reasonable for all of the solutes studied here, since the number of hydrogen-bonded water molecules only begins to drop at solute concentrations above 3M (Figs. 1 and 4).

Another important assumption is that hydration water moves with the solute, while bulk water moves more like the pure liquid. In the case of TMAO, the bulk moves significantly more slowly than the pure liquid but not in concert with the solute i, and its diffusion rate can be predicted from the concentration dependence of Di(c). To examine how water is behaving in the simulation, DhydMD(c), the average diffusion coefficient for water that is hydrogen bonded to a solute, and DbulkMD(c), the average diffusion coefficient for water that is not hydrogen bonded to a solute (i.e., bulk water), are also calculated from the simulations (see Sec. II).

First, the diffusion of water hydrogen bonded to the solute is examined. When the bulk water diffusion can be assumed to be independent of solute concentration (i.e., ai = 0), then DhydMD(c) is predicted well by Di(0). For instance, DhydMD(1M) is 1.24 × 10−9 m2/s and Di(0) is 1.18 × 10−9 m2/s in NaCl solution and DhydMD(1M) is 1.53 × 10−9 m2/s and Di(0) is 1.59 × 10−9 m2/s in KCl solution. However, when the bulk water diffusion is strongly dependent on the solute concentration, DhydMD(c) is predicted better by Di(c) which becomes more important at higher concentrations. In TMAO solution, DhydMD(3M) is 0.43 × 10−9 m2/s, while Di(0) is 0.94 × 10−9 m2/s and Di(c) is 0.38 × 10−9 m2/s, so it is better predicted by Di(c). In urea solution, DhydMD(3M) is 1.47 × 10−9 m2/s, while Di(0) is 1.73 × 10−9 m2/s and Di(c) is 1.47 × 10−9 m2/s × 10−9 m2/s, so it is also better predicted by Di(c), but the difference between Di(0) and Di(c) is relatively small and it is difficult to get an accurate value of ai when the slope is small.

Next, the dependence of the bulk water diffusion on solute concentration is examined. For the most part, DbulkMD(c) is predicted well by Dw(0) at low concentrations. For instance, DbulkMD(1M) is 2.22 × 10−9 m2/s in NaCl solution and 2.37 × 10−9 m2/s in KCl solution, which are both comparable with the pure liquid value of Dw = 2.62 × 10−9 m2/s for TIP4P-Ew in a similar size box. In addition, since DbulkMD(3M) is 2.21 × 10−9 m2/s in urea solution, it compares well with Dw = 2.80 × 10−9 m2/s, the pure liquid value for TIP4P-Ew in a similar size box. However, since DbulkMD(3M) is 1.28 × 10−9 m2/s in TMAO solution, it is better predicted when a non-zero value is used, which gives Dbulk(c) = 1.38 × 10−9 m2/s.

Overall, in the systems studied here, the concentration dependence of Dbulk(c) is largest in the TMAO solutions. In particular, the marked decrease in DTMAO(c) with concentration appears to be due to the significantly increasing viscosity of the bulk liquid, indicating that Dbulk(c) decreases with concentration. On the other hand, the concentration dependence of Dbulk(c) is smallest of the systems studied here in the urea solutions, indicating that Dbulk(c) stays close to Dw(0), the value for pure water. This indicates that TMAO may act as a kosmotrope as far as it has effects on protein stability by increasing the viscosity of water. Adding zwitterionic and polar solutes might increase the viscosity due to dipole ordering (the dipole moments of TMAO and urea are 6.68 D and 5.11 D, respectively, using CHARMM parameters), which leads to large Kirkwood g-factors.45 Although not tested, the greater effect of TMAO may be because the dipole is largely linear, while the positive charge in urea is spread between the two nitrogens. Interestingly, since the HB-B model for the ternary solution uses hydrogen bond numbers and diffusion coefficients for TMAO and urea obtained from the binary solutions rather than those from the ternary solution simulations, the good agreement for the ternary mixtures indicates that effects of TMAO and urea on bulk water diffusivity are additive.

Further examination of the salt solution results indicates that the small deviations between the calculated or experimental Dw(c) and the HB-B results may lie in the bulk water. For instance, the simulation results for NaCl solutions show that the HB-B model predicts a smaller decrease in Dw(c) with concentration than calculated from the simulation [Fig. 3(a)]. However, Na+ might be expected to slow down water beyond the first shell since the waters directly hydrogen bonded to the ion have a long lifetime of ∼37 ps and thus may also slow down other waters hydrogen bonded to them. Water molecules that are not directly bonded to the Na+ are counted as “bulk,” so the decrease in Dbulk(c) would be underpredicted. In addition, the experimental results for KCl indicate a slight initial increase in Dw(c) with concentration, which might indicate a decrease in the water structure so that the bulk would diffuse faster with increasing solute concentrations, which is not seen in the HB-B model (Fig. 3). Overall, the differences in agreement between Dw(c) and the HB-B model Dw(c) (Fig. 3) between simulations and experiments imply that bulk TIP4P-Ew water might have differences in behavior from real bulk water. In addition, comparison of Dw(c) from the HB-B model and from direct calculations for the simulation suggests that another important factor in predicting the anomalous diffusion of water in electrolyte solutions as a function of ionic concentration might be the number and lifetime of water molecules directly hydrogen bonded to the ions, which affects the average diffusion of water directly by the number diffusing at the ion rate and indirectly by affecting the diffusion rate of the ions. In particular, according to the HB-B model, the too slow rate of ion diffusion in the simulations gives rise to the too steep decrease in Dw(c) with ionic concentration compared to experiment since the number of hydrogen bonded water molecules was essentially the same in simulations and experiments (see Table S6 in the supplementary material).

Thus, the HB-B model appears to be a good description of the effects of ions, zwitterions, and polar molecules on the diffusion of water in aqueous solutions in terms of diffusion of hydration water and of bulk water. The results so far indicate that the importance of accounting for the diffusion of hydrogen-bonded water appears most important for ions, while the diffusion of bulk water is most important for TMAO possibly due to long range ordering. In addition, since the water diffusion in salt solutions is predicted well by low concentration values, the HB-B model can be used to estimate the high concentration behavior for optimization of force fields.

In its simplest form, the HB-B model (with ai = 0) describes Dw(c) at high concentrations in terms of quantities that can be obtained at very low concentration, namely, niw(0) and Di(0) as well as Dw(0).

The diffusivities in aqueous solutions of NaCl, KCl, TMAO, urea, and TMAO-urea are examined in molecular dynamics simulations. The decrease in the diffusivity of water with concentration can be described by the HB-B model in which the water molecules hydrogen bonded to the solutes are assumed to diffuse at the same rate as the solutes, while the remainder of the water molecules diffuse at a rate corresponding to a bulk liquid. In the simulations of the ionic and urea solutions, the diffusivity of bulk water can be described as that of pure water so that the major factor determining the decrease in diffusion with concentration appears to be the number of waters hydrogen bonded to the ion. On the other hand, in the simulations of TMAO solutions, the diffusivity of bulk water apparently decreases with increasing solute concentrations, presumably because TMAO produces more order in water beyond the directly hydrogen-bonded water. This change in the diffusivity of bulk water appears to be more important in determining the diffusivity of all water. Although the specific results here may be force-field dependent, the HB-B model provides a useful framework for understanding the effects of hydration shell water versus bulk water, whether in simulations or in experiments. In particular, inconsistencies between simulation and experiment may help identify deficiencies in force fields.

See supplementary material for data for the simulation (Tables S1 to S5) and the HB-B model (Tables S6 to S7). Plots are also given for the coordination numbers of water around ions as a function of ion concentration (Fig. S1) as well as for the number of hydrogen bonds (Fig. S2) and solute diffusion coefficients (Fig. S3) in TMAO and urea solutions as a function of total solute concentration.

T.I., C.C.D., and Q.H. acknowledge support from the National Science Foundation through Grant No. CHE-1464766, X.T. acknowledges support from the National Institutes of Health through Grant No. R01-GM122441, and Q.H. acknowledges support from the National Institutes of Health through Grant No. R21-GM104500 and from the McGowan Charitable Fund. This work used time on the LoBoS cluster at the Laboratory for Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, which was generously provided by Dr. Bernard R. Brooks, the Extreme Science and Engineering Discovery Environment (XSEDE) granted via No. MCB990010, which is supported by National Science Foundation Grant No. OCI-1053575, and the Medusa cluster maintained by University Information Services at Georgetown University.

1.
K. J.
Müller
and
H. G.
Hertz
,
J. Phys. Chem.
100
,
1256
(
1996
).
2.
P. B.
Ishai
,
E.
Mamontov
,
J. D.
Nickels
, and
A. P.
Sokolov
,
J. Phys. Chem. B
117
,
7724
(
2013
).
3.
R.
Mills
and
V. M. M.
Lobo
,
Self-Diffusion in Electrolyte Solutions
(
Elsevier
,
Amsterdam
,
1989
).
4.
J. S.
Kim
,
Z.
Wu
,
A. R.
Morrow
,
A.
Yethiraj
, and
A.
Yethiraj
,
J. Phys. Chem. B
116
,
12007
(
2012
).
5.
Y.
Ding
,
A. A.
Hassanali
, and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
111
,
3310
(
2014
).
6.
Y.
Yao
,
Y.
Kanai
, and
M. L.
Berkowitz
,
J. Phys. Chem. Lett.
5
,
2711
(
2014
).
7.
Y.
Yao
,
M. L.
Berkowitz
, and
Y.
Kanai
,
J. Chem. Phys.
143
,
241101
(
2015
).
8.
A. J.
Lee
and
S. W.
Rick
,
J. Chem. Phys.
134
,
184507
(
2011
).
9.
M.
Soniat
and
S. W.
Rick
,
J. Chem. Phys.
137
,
044511
(
2012
).
10.
Z. R.
Kann
and
J. L.
Skinner
,
J. Chem. Phys.
141
,
104507
(
2014
).
11.
Y.
Marcus
,
Chem. Rev.
109
,
1346
(
2009
).
12.
D. J.
Tobias
and
J. C.
Hemminger
,
Science
319
,
1197
(
2009
).
13.
E. J.
Guinn
,
L. M.
Pegram
,
M. W.
Capp
,
M. N.
Pollock
, and
M. T.
Record
, Jr.
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
16932
(
2011
).
14.
L.
Larini
and
J.-E.
Shea
,
J. Phys. Chem. B
117
,
13268
(
2013
).
15.
C.
Hölzl
,
P.
Kibies
,
S.
Imoto
,
R.
Frach
,
S.
Suladze
,
R.
Winter
,
D.
Marx
,
D.
Horinek
, and
S. M.
Kast
,
J. Chem. Phys.
144
,
144104
(
2016
).
16.
S.
Weerasinghe
and
P. E.
Smith
,
J. Phys. Chem. B
107
,
3891
(
2003
).
17.
L. J.
Smith
,
H. J. C.
Berendsen
, and
W. F.
van Gunsteren
,
J. Phys. Chem. B
108
,
1065
(
2004
).
18.
J. A.
Schellman
,
Biophys. Chem.
96
,
91
(
2002
).
19.
A.
Idrissi
,
Spectrochim. Acta, Part A
61
,
1
(
2005
).
20.
M. C.
Stumpe
and
H.
Grubmüller
,
J. Phys. Chem. B
111
,
6220
(
2007
).
21.
P. H.
Yancey
,
M. E.
Clark
,
S. C.
Hand
,
R. D.
Bowlus
, and
G. N.
Somero
,
Science
217
,
1214
(
1982
).
22.
T. Y.
Lin
and
S. N.
Timasheff
,
Biochemistry
33
,
12695
(
1994
).
23.
P. H.
Yancey
,
A. L.
Fyfe-Johnson
,
R. H.
Kelly
,
V. P.
Walker
, and
M. T.
Auñón
,
J. Exp. Zool.
289
,
172
(
2001
).
24.
P. H.
Yancey
,
M.
Gerringer
,
A. A.
Rowden
,
J. C.
Drazen
, and
A.
Jamieson
,
Proc. Natl. Acad. Sci. U. S. A.
111
,
4461
(
2014
).
25.
Q.
Zou
,
B. J.
Bennion
,
V.
Daggett
, and
K. P.
Murphy
,
J. Am. Chem. Soc.
124
,
1192
(
2002
).
26.
S.
Paul
and
G. N.
Patey
,
J. Am. Chem. Soc.
129
,
4476
(
2007
).
27.
B. R.
Brooks
,
C. L.
Brooks
 III
,
A. D.
MacKerell
, Jr.
,
L.
Nilsson
,
R. J.
Petrella
,
B.
Roux
,
Y.
Won
,
G.
Archontis
,
C.
Bartels
,
S.
Boresch
,
A.
Caflisch
,
L.
Caves
,
Q.
Cui
,
A. R.
Dinner
,
M.
Feig
,
S.
Fischer
,
J.
Gao
,
M.
Hodoscek
,
W.
Im
,
K.
Kuczera
,
T.
Lazaridis
,
J.
Ma
,
V.
Ovchinnikov
,
E.
Paci
,
R. W.
Pastor
,
C. B.
Post
,
J. Z.
Pu
,
M.
Schaefer
,
B.
Tidor
,
R. M.
Venable
,
H. L.
Woodcock
,
X.
Wu
,
W.
Yang
,
D. M.
York
, and
M.
Karplus
,
J. Comput. Chem.
30
,
1545
(
2009
).
28.
B. T.
Miller
,
R. P.
Singh
,
J. B.
Klauda
,
M.
Hodošček
,
B. R.
Brooks
, and
H. L.
Woodcock
 III
,
J. Chem. Inf. Model.
48
,
1920
(
2008
).
29.
H. W.
Horn
,
W. C.
Swope
,
J. W.
Pitera
,
J. D.
Madura
,
T. J.
Dick
,
G. L.
Hura
, and
T.
Head-Gordon
,
J. Chem. Phys.
120
,
9665
(
2004
).
30.
I. S.
Joung
and
T. E.
Cheatham
 III
,
J. Phys. Chem. B
112
,
9020
(
2008
).
31.
A. D.
MacKerell
, Jr.
,
D.
Bashford
,
M.
Bellot
,
R. L.
Dunbrack
, Jr.
,
M. J.
Field
,
S.
Fischer
,
J.
Gao
,
H.
Guo
,
S.
Ha
,
D.
Joseph
,
K.
Kuchnir
,
K.
Kuczera
,
F. T. K.
Lau
,
M.
Mattos
,
S.
Michnick
,
D. T.
Nguyen
,
T.
Ngo
,
B.
Prodhom
,
B.
Roux
,
M.
Schlenkrich
,
J.
Smith
,
R.
Stote
,
J.
Straub
,
J.
Wiorkiewicz-Kuczera
, and
M.
Karplus
,
J. Phys. Chem. B
102
,
3586
(
1998
).
32.
J. P.
Rychaert
,
G.
Ciccotti
, and
H. J. C.
Berendsen
,
J. Comput. Phys.
23
,
327
(
1977
).
33.
L.
Verlet
,
Phys. Rev.
165
,
201
(
1968
).
34.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
35.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
36.
K. N.
Tran
,
M.-L.
Tan
, and
T.
Ichiye
,
J. Chem. Phys.
145
,
034501
(
2016
).
37.
I.-C.
Yeh
and
G.
Hummer
,
J. Phys. Chem. B
108
,
15873
(
2004
).
38.
J.
Kestin
,
H. E.
Khalifa
, and
R. J.
Correia
,
J. Phys. Chem. Ref. Data
10
,
71
(
1981
).
39.
J.
Kestin
,
H. E.
Khalifa
, and
R. J.
Correia
,
J. Phys. Chem. Ref. Data
10
,
57
(
1981
).
40.
J.
Hunger
,
N.
Ottosson
,
K.
Mazur
,
M.
Bonn
, and
H. J.
Bakker
,
Phys. Chem. Chem. Phys.
17
,
298
(
2015
).
41.
K.
Kawahara
and
C.
Tanford
,
J. Biol. Chem.
241
,
3228
(
1966
).
42.
R.
Mancinelli
,
A.
Botti
,
F.
Bruni
,
M. A.
Ricci
, and
A. K.
Soper
,
J. Phys. Chem. B
111
,
13570
(
2007
).
43.
R.
Sinbaldi
,
C.
Casieri
,
S.
Mechioma
,
G.
Onori
,
A. L.
Segre
,
S.
Viel
,
L.
Mannina
, and
F.
De Luca
,
J. Phys. Chem. B
110
,
8885
(
2006
).
44.
M.
Mayele
and
M.
Holz
,
Phys. Chem. Chem. Phys.
2
,
2429
(
2000
).
45.
J. G.
Kirkwood
,
J. Chem. Phys.
7
,
911
(
1939
).

Supplementary Material