Hydration and self-aggregation of a neutral cosolute from dielectric relaxation spectroscopy and MD simulations: the case of 1,3-dimethylurea

Abstract

The influence of the amphiphile 1,3-dimethylurea (1,3-DMU) on the dynamic properties of water was studied using dielectric relaxation spectroscopy. The experiment provided evidence for substantial retardation of water reorientation in the hydration shell of 1,3-DMU, leading to a separate slow-water relaxation in addition to contributions from bulk-like and fast water as well as from the solute. From the amplitudes of the resolved water modes effective hydration numbers were calculated, showing that each 1,3-DMU molecule effectively freezes the reorientation of 1–2 water molecules. Additionally, a significant amount of solvent molecules, decreasing from ∼39 at infinite dilution to ∼3 close to the solubility limit, is retarded by a factor of ∼1.4 to 2.3, depending on concentration. The marked increase of the solute amplitude indicates pronounced parallel dipole alignment between 1,3-DMU and its strongly bound H₂O molecules. Molecular dynamics (MD) simulations of selected solutions revealed a notable slowdown of water rotation for those solvent molecules surrounding the methyl groups of 1,3-DMU and strong binding of ∼2H₂O by the hydrophilic carbonyl group, corroborating thus the experimental results. Additionally, the simulations revealed 1,3-DMU self-aggregates of substantial lifetime.

---

Introduction

Various aspects of hydration in aqueous solutions of urea and its derivatives are of great interest, particularly when related to their denaturing properties. The denaturation of proteins by urea itself has been controversial for a long time, and despite decades of research¹ agreement on a direct mechanism, driven by preferential binding of urea, has emerged only in recent years.^2–6 Alkyl and aryl substituted ureas did not receive this attention despite their applications in biochemistry as even stronger denaturants than urea,^2,7,8 and while an indirect effect – i.e. via modulation of the water structure – has been ruled out in the case of the parent compound, its possible contribution to the denaturation power of N-substituted urea derivatives has not been studied.

Similar to urea, 1,3-dimethylurea (1,3-DMU) is used for the denaturation of proteins and the investigation of their structure.^9–11 Like its parent, this compound has carbonyl and amino groups which can directly hydrogen bond to hydrophilic parts of proteins.² Furthermore, urea can also directly bind to nonpolar groups by van der Waals interactions,¹² and the presence of alkyl groups will render this interaction even more favourable.

Indirect cosolvent effects on protein stability are intricately linked to the thermodynamic properties of the cosolvent–water mixture and thus originate from the underlying molecular-level structure. For urea, indirect effects have been excluded, because the molecule fits nicely into the hydrogen bond network of water.^13,14 This is unlikely for 1,3-DMU molecules because of the presence of two methyl groups which are expected to significantly affect the static and dynamic properties of the solvent. Indeed, a significant slowdown of water dynamics around the methyl groups of various osmolytes, such as tetramethylurea (TMU),^15–17 trimethylamine-N-oxide (TMAO),^16,18 proline,¹⁶ and N-methylacetamide,¹⁶ was observed. Accordingly, some authors^7,19–21 claimed indirect contributions to be important in protein denaturation by substituted ureas.

Recently Sagle et al.²² investigated the conformational changes of the thermo-responsive polymer poly(N-isopropylacrylamide) (PNIPAM) in the presence of urea and its methyl substituted derivatives. No evidence for direct binding of 1,3-DMU to the amide moieties of the polymer was found. This contrasts the behavior of urea which directly interacts with the polymer by cross-linking the amide groups via hydrogen bonding. Molecular dynamics (MD) simulations²³ of aqueous PNIPAM solutions, both in the presence and the absence of urea or 1,3-DMU, evidenced notable 1,3-DMU self-aggregation which apparently favours the unfolded state of PNIPAM and increases the lower critical solution temperature (LCST) of the mixture. On the other hand, urea binding to the folded state of PNIPAM decreased the LCST.

In view of the above findings, it is of interest to understand how 1,3-DMU behaves in aqueous solution. In a few experimental studies it was shown that this solute has a distinct effect on the dynamics of water molecules in its hydration shell. In particular, NMR experiments²⁴ revealed a downfield shift of the water proton signal in the presence of 1,3-DMU, indicating an enhancement of the structure of the hydrating water molecules. Dielectric relaxation spectroscopy (DRS)²⁵ indicated retarded rotational dynamics of hydration water due to the presence of hydrophobic methyl groups. The idea of a strengthened water structure around 1,3-DMU is also supported by data for heat capacities,^24,26 molar volumes,^26,27 and molar expansion coefficients.²⁷ Nevertheless, these findings give no information on the magnitude and the molecular origin of the retarded hydration-shell dynamics.

As mentioned above, a comprehensive picture of their hydration pattern in aqueous solutions is a prerequisite for understanding the denaturation properties of substituted ureas. Among the experimental techniques capable of probing the hydration of polar molecular solutes, dielectric relaxation spectroscopy yields an excitingly wide scope of information both on structure and dynamics.²⁸ Recently, this technique was applied to study aqueous solutions of urea²⁹ and cosolutes such as TMU,¹⁵ TMAO,^18,30 and ectoine.³¹ While urea strongly binds ∼2H₂O molecules (where the term “strong binding” refers to the retardation of solvent molecules to such an extent that their contribution coincides with the solute) but otherwise has only small effects on water dynamics, it was possible for the other compounds to distinguish and quantify two fractions of bound water with significantly different dynamics. The disadvantage of DRS here is that no direct, i.e. molecular-level, information on the solute interaction sites of these two water fractions can be obtained. In this regard, molecular dynamics (MD) simulations provide valuable complementary insights in addition to a broad range of other important data. Accordingly, this contribution summarizes the results of a joint DRS and MD study on aqueous solutions of 1,3-DMU at 25 °C, covering solute molalities up to m = 8.0 mol kg⁻¹. For the simulations an especially derived force-field model for the solute was used.

Method section

Experimental

Dielectric relaxation spectroscopy probes the macroscopic polarization of the sample in a time-dependent electric field.³² In the case of solutions of non-electrolytes, polarization originates from the reorientation of permanent dipoles (both solvent and solute molecules) and from intramolecular polarizability. As a function of frequency, ν, the response can be expressed in terms of the complex permittivity

[small epsilon, Greek, circumflex](ν) = ε′(ν) − iε′′(ν), (1)

where the relative permittivity, ε′(ν), expresses how far dipoles can follow the electric field and the dielectric loss, ε′′(ν), is a measure of the energy dissipated in the sample. With increasing frequency ε′(ν) drops from the static relative permittivity, ε = ε′(ν → 0), to the high-frequency limit, ε_∞, characteristic of intramolecular polarizability.

The present dielectric spectra were recorded at (25.00 ± 0.05) °C in the frequency range of 0.2 ≤ ν/GHz ≤ 89. At ν ≤ 50 GHz a frequency-domain reflectometer, based on Agilent 85070E-20 and 85070E-50 dielectric probes connected to an Agilent E8364B vector network analyzer (VNA) and an electronic calibration module (Agilent N4693B), was used.³³ A range of 60–89 GHz was covered using a waveguide interferometer.³⁴ For selected samples the reliability of the reflection measurements was additionally crosschecked with a 27–40 GHz variable-pathlength waveguide transmission cell hooked to the VNA. Examples of the spectra are shown in Fig. 1 and 2.

Fig. 1 (A) Relative permittivity, ε′, and (B) dielectric loss, ε′′, spectra of aqueous 1,3-DMU solutions at 25 °C and c/M = 0.1157; 0.2519; 0.4794; 0.9156; 1.7079; 2.4521; 2.9910; 4.1322; 4.9118. Symbols show experimental data (omitted for some samples for clarity); lines show the fit using eqn (2); arrows indicate the trend of increasing solute concentration.

Fig. 2 Relative permittivity, ε′(ν), and dielectric loss, ε′′(ν), spectrum of 3.0 M aqueous 1,3-DMU at 25 °C. The symbols indicate experimental values; lines show the results of the fit with eqn (2), with shaded areas indicating the loss of contributions of the resolved modes.

Aqueous solutions of 1,3-DMU were prepared by weight without buoyancy correction from degassed Millipore Milli-Q water and 1,3-DMU (Sigma-Aldrich, >99%) vacuum-dried at 50 °C to constant mass. Densities, d, necessary for calculating molar solute concentrations, c [in M ≡ mol L⁻¹], were determined using an Anton Paar DMA 5000 M vibrating tube densimeter yielding d with a nominal uncertainty of ±5 × 10⁻⁶ g cm⁻³ at (25.000 ± 0.005) °C. Solution viscosity, η, was measured at (25.00 ± 0.01) °C with an Anton Paar AMVn rolling ball viscometer using a 1.6 mm capillary calibrated with degassed Millipore Milli-Q water. The densities and viscosities of the solutions are given in Table S1 of the ESI.†

The dielectric spectra were evaluated by simultaneously fitting ε′(ν) and ε′′(ν) to conceivable relaxation models based on sums of n, n = 1…5, Havriliak–Negami equations or their simplified variants and the resulting fits and their parameters were scrutinized as described in detail elsewhere.^35,36 It was found that all solution spectra were best fitted by the sum of four Debye equations (the 4D model), i.e.

where – from low to high frequencies – S_j and τ_j are the relaxation amplitudes and relaxation times of the resolved processes for 1,3-DMU, slow, bulk-like and fast water; see below for the assignment. The obtained parameters are summarized in Table S2 of the ESI;† examples of the fits are shown in Fig. 1 and 2.

Molecular dynamics simulations

The equilibrium geometry of an isolated 1,3-DMU molecule was determined from quantum-chemical calculations utilizing the Gaussian 03 program.³⁷ Calculations were performed at the density functional theory (DFT) level using the B3LYP exchange–correlation functional³⁸ combined with the cc-PVTZ basis set.³⁹ Implicit solvent-medium effects on the ground state geometry were taken into account by the self-consistent reaction field method (SCRF) via the self-consistent isodensity polarizable continuum model (PCM)⁴⁰ implemented in the software. Default parameter values were used. The atomic charge distribution was obtained using the RESP-fit method.^41,42 The obtained partial charges and structural characteristics of 1,3-DMU, together with the corresponding force-field parameters of all atoms of the molecule, namely carbonyl carbon, C, oxygen, O, nitrogen, N, hydrogen, H_N, methyl-group carbon, C_H, and methyl-group hydrogen, H_C, are given in the ESI.†

The equilibrium values for the bond lengths, angles and torsions used in the simulations were taken from the above DFT calculations, whereas the used force constants are those of the GROMOS 53 A6 force field version.⁴³ Non-bonding van der Waals interactions were calculated using the Lennard-Jones 12-6 potential. The homoatomic pair parameters of the KBFF urea model^44,45 were employed, whereas heteroatomic pair parameters were calculated by the geometric combination rule. Lacking Lennard-Jones parameters for the methyl carbons, C_H, were calculated using the procedure given in ref. 45 based on a correlation between the atomic size and the atomic hybrid components of molecular polarizabilities.⁴⁶ Required values for the atomic hybrid components and ionization potentials were taken from ref. 46 and 47, respectively. For the methyl hydrogen atoms, H_C, the Lennard-Jones parameters were obtained by fitting the experimental values of the density and the diffusion coefficient of a 3.0 M solution. The obtained density and self-diffusivity values given in Table 1 show good agreement with experiment.^48–50

Table 1 Selected properties of the simulated systems: numbers of 1,3-DMU, N₁, and water molecules, N₂, in the simulation box, corresponding to solute molality, m/mol kg⁻¹, molarity, c/M, simulated (calc.) and experimental (exp.) density, d/g cm⁻³, and solute diffusion coefficient, D₁ × 10⁵/cm² s⁻¹

N 1	N 2	m	c	d		D 1
				Calc.	Exp.^b	Calc.	Exp.
a Obtained using the density values, dcalc, found using the present MD simulations. b Data of ref. 48 interpolated with a third-order polynomial. c Data of ref. 49 interpolated with polynomials. d Data of ref. 50 interpolated with polynomials.
26	2921	0.4941	0.4755	1.0043	1.0012	1.3 ± 0.5	0.9144^c; 0.8630^d
51	2814	1.0060	0.9327	1.0093	1.0053	1.0 ± 0.1	0.8028^c; 0.7757^d
95	2625	2.0089	1.7373	1.0179	1.0130	0.89 ± 0.08	0.6750^c; 0.6631^d
167	2320	3.9957	3.0463	1.0308	1.0261	0.56 ± 0.02	0.5828^c; 0.5342^d
225	2085	5.9902	4.0780	1.0401	1.0367	0.5 ± 0.1	0.5352^c; 0.4706^d
274	1899	8.0092	4.9172	1.0472	1.0447	0.42 ± 0.01	0.4629^c; 0.4493^d

---

MD simulations were performed by means of the GROMACS software package (version 4.5.6).⁵¹ An isothermal–isobaric (NPT) ensemble with periodic boundary conditions and an integration time step of 1 fs was used. The temperature was set to 25 °C and controlled using a stochastic v-rescaling thermostat⁵² using a time constant of 0.1 ps for heat bath coupling. The pressure was controlled using an isotropic Parrinello–Rahman barostat^53,54 at 1 bar using a time constant of 1 ps for pressure relaxation. The van der Waals interactions were computed with a cutoff of 1.0 nm with a smooth decay of the forces starting at 0.9 nm. Electrostatic interactions were computed by means of the particle-mesh Ewald technique^55,56 with maximum space grid points being 0.12 nm.

Cubic boxes with lengths of 4.5 nm were created with water and solute molecules evenly distributed in the starting configuration and systems were then equilibrated for 5 ns. Depending on the desired properties, two different runs were performed: (i) structural properties were calculated from the configurations generated during the subsequent 45 ns, with coordinates being saved every 5000 steps; (ii) to calculate contact and H-bond dynamics, a 2 ns run was used where the coordinates were saved every 200 steps. Analysis was conducted using the GROMACS and TRAVIS⁵⁷ programs. A summary of the performed simulations is given in Table 1.

Results and discussion

Dielectric relaxation spectroscopy

Representation of dielectric spectra. Although the solution spectra showed only a single loss peak, albeit of strongly varying shape (Fig. 1), all of them were best fitted by the 4D model, eqn (2). An example for this decomposition, indicating the individual contributions, is shown in Fig. 2. However, in free-running fits the amplitude of the lowest-frequency mode, S₁, scattered considerably for c ≤ 0.5 M due to its small magnitude and location close to the low-frequency limit of the experiment. Therefore, its relaxation time, τ₁, was fixed to values extrapolated from a straight-line fit to the corresponding data at higher concentrations (Table S2 of the ESI†).

The two higher-frequency modes, (S₃, τ₃) and (S₄, τ₄), are readily assigned to the solvent as they also appear for pure water⁵⁸ and – shown in Fig. 3 – their parameters smoothly evolve from the values at c = 0, obtained in our laboratory from fitting [small epsilon, Greek, circumflex](ν) in the 0.1–2000 GHz range (S₃(0) = 72.42, τ₃(0) = 8.35 ps, S₄(0) = 2.43, τ₄(0) = 0.278 ps, ε_∞(0) = 3.52).⁵⁹ The first peak, centred at ∼18 GHz, is due to the cooperative relaxation of the hydrogen-bond network of bulk-like water (i.e. more or less unperturbed H₂O molecules behaving similar to the neat solvent). The underlying mechanism can be envisaged along the lines of the Laage–Hynes model of water dynamics⁶⁰ as a wait, jump and settle sequence for the water dipole and its surroundings. According to this picture τ₃ is essentially determined by the waiting time for the water molecule between two fast large-angle jumps, with this waiting time strongly influenced by the time required by the environment to equilibrate with the new dipole orientation. The fast dipole jump of this mechanism, leading to a change of H-bond allegations, gives rise to small-amplitude second (fast) water relaxation, at ∼500 GHz. The lowest-frequency mode (S₁, τ₁), shifting from ∼4 to ∼1.5 GHz with increasing concentration, can be attributed to the solute, whereas the relaxation at ∼10 GHz, the “slow-water” mode (S₂, τ₂), is typical for retarded H₂O molecules in contact with hydrophobic moieties^18,61 but appears also with solutes forming H-bonds somewhat stronger than those among water molecules.^59,62 While S₁ exhibits a linear increase with increasing solute molarity, S₂ passes through a maximum at c ≈ 1.7 M (Fig. 3). The retardation factor for this slow water fraction increases from r = τ₂/τ₃ = 1.4 at the lowest solute concentration to 2.3 at c = 4.9 M.

Fig. 3 (A) Amplitudes, S_j, and (B) relaxation times, τ_j, of the 1,3-DMU (j = 1), slow (j = 2), bulk (j = 3), and fast water (j = 4) modes obtained from a fit of the spectra with the 4D model. Symbols indicate the values of Table S2 (ESI†); lines are low-order polynomial fits.

A remarkable feature of the present spectra, irrespective of the chosen formal description and in line with previous findings,²⁵ is the monotonic increase of ε′ at ν < 0.5 GHz with increasing solute concentration (Fig. 1A). Thus, also the static permittivity, (Table S2 of the ESI†), increases. This observation is surprising because, based on the respective analytical concentrations, a decrease of ε from the pure-water value of 78.4 to 61.1 at c = 4.9 M would be expected from the effective dipole moments, μ_eff, of 1,3-DMU (5.5 D; from quantum-chemical calculations) and water (3.84 D⁶³). This indicates that the dielectric relaxation of 1,3-DMU–water mixtures does not simply arise from independent reorientations of solute and solvent dipoles. The reasons for that will be explored in the following sections.

The starting point for the evaluation of the obtained relaxation amplitudes, S_j, and also of the above estimation of ε is the equation^59,64

which relates S_j to the molar concentration, c_j, and the effective dipole moment, μ_j,eff, of species j. For 1,3-DMU a cavity-field factor of A_j = 0.223 was estimated from its geometric parameters, whereas A_j = 1/3 was taken for the water modes. Note that μ_j,eff = g_j^1/2μ_j,ap where g_j is the equivalent to the Kirkwood factor accounting for orientational correlations among the dipoles. The apparent dipole moment in solution, μ_j,ap = μ_j/(1 − f_jα_j), arises from the permanent (gas-phase) dipole moment, μ_j, of the molecule and a reaction-field correction (field factor f_j) due to its polarizability, α_j.⁶⁴

Water relaxation and hydration numbers. With relaxation times 0.46 ≤ τ₄/ps ≤ 3.12 (Table S2, ESI†) the fast water mode is essentially outside the current frequency range and just its low-frequency wing is detected. Accordingly, the obtained values of S₄ and τ₄ are rather noisy and we will refrain from discussing a possible concentration dependence of the relaxation time. On the other hand, the relaxation time of the cooperative mode of bulk-like water, τ₃, is well defined and its value slightly increases with c (Fig. 3B; Table S2 of the ESI†). It is also roughly proportional to solution viscosity, η (Fig. S1 of the ESI†), but with a rather small slope, dτ₃/dη = (0.74 ± 0.14) 10⁻¹¹ Pa⁻¹, which excludes – as expected⁶⁰ – the rotational diffusion of individual water molecules as the relaxation mechanism. Initially, the corresponding amplitude, S₃, sharply drops with increasing concentration but remains rather constant at c ≥ 2.5 M (Fig. 3A).

According to the Laage–Hynes model of water dynamics,⁶⁰ and in line with previous observations,⁶¹ the bulk-like mode S₃ and fast mode S₄ do not represent different water species with distinct dynamics but merely highlight the two dominating steps in the time evolution of water reorientation. For this reason both amplitudes were combined with the amplitude of “unperturbed” (not bound to solute) water

S_un(c) = S₃(c) + S₄(c) + ε_∞(c) − ε_∞(0) (4)

for further analysis. Note that the last two terms of eqn (4) correct to some extent for the lack of terahertz data in the solution spectra, which results in somewhat too large and noisy ε_∞(c) values.

With the help of eqn (3), normalized to pure water,⁶¹ the concentrations of water not bound to solute, c_un, and slow water, c_s, were then calculated from S_un(c) and S₂(c), respectively. In conjunction with the analytical water concentration, c_w, the obtained results then yielded the effective hydration numbers

Z_t = (c_w − c_un)/c, (5)

Z_s = c_s/c (6)

and

Z_ib = Z_t − Z_s, (7)

where Z_t is the total number of H₂O molecules per equivalent of solute differing in their dynamics from bulk-like water; Z_s is the number of retarded (slow) but still DRS-detected solvent molecules and Z_ib represents the irrotationally bound (“frozen”, not DRS-detected) water.⁵⁹Fig. 4 shows the obtained values for Z_t, Z_s and Z_ib as a function of solute concentration.

Fig. 4 Numbers of frozen, Z_ib, slow, Z_s, and total, Z_t = Z_ib + Z_s, water molecules per 1,3-DMU from DRS and number of H₂O molecules surrounding the solute, CN_MD(H₂O), obtained with MD simulations (see the details below) as a function of 1,3-DMU molarity at 25 °C. Symbols represent experimental data; lines are appropriate fits (open symbols not included); the black dashed line indicates Z_i = 0; the vertical magenta dashed line corresponds to the molarity at which the first hydration shells of 1,3-DMU start to overlap.

As generally found for hydrophobic solutes,^18,31,61,65 the present slow-water hydration number, Z_s, strongly decreases with increasing solute concentration, dropping from ∼39 at c → 0 to ∼3 at 4.9 M. This is almost certainly due to the hydration-shell overlap, because at the highest solute concentration the solute : solvent molar ratio has dropped to 1 : 7. Interestingly the dilute-solution values of Z_s for 1,3-DMU are significantly larger than those for TMU¹⁵ (∼12), TMAO¹⁸ (∼9) and ectoine³¹ (∼12). On the other hand, the amount of strongly bound water, Z_ib ≈ 1.6, is independent of solute concentration (Fig. 4) and significantly smaller than the number of H₂O molecules in the first hydration shell determined by MD which drops from ∼25.9 to 12.5⁶⁶ (Fig. 4).

Solute relaxation. In line with previous studies²⁵ the lowest-frequency relaxation, peaking at (∼1.5 to 4) GHz depending on c (Fig. 3), was assigned to the solute. However, there are strong arguments that this mode is not solely due to the reorientation of bare 1,3-DMU dipoles.

One argument comes from the length of its relaxation time, τ₁, which increases from 37 ps at c = 0.12 M to 115 ps at the highest concentration of this study, 4.9 M. The corresponding single-molecule rotational correlation times, τ_rot,1, calculated from τ₁ using the Powles–Glarum equation^67,68

where ε₁ = ε − S₁, yielded a straight line of slope a = (2.78 ± 0.06) (Pa s)⁻¹ when plotted against solution viscosity, η (Fig. S2 of the ESI†), i.e. followed the extended Stokes–Einstein–Debye equation.⁶⁹ This finding suggests the rotational diffusion of the causing dipolar entity with an effective volume of V_eff = ak_BT/3 = (38.1 ± 0.8) ų. At the same time, quantum chemical calculations yielded a molecular volume of V_m = 122 ų for 1,3-DMU. According to Dote et al.,⁶⁹V_eff is related to V_m through the expression V_eff = f_⊥CV_m, where the factor f_⊥ (=1.23 for 1,3-DMU) accounts for the deviation of molecular geometry from the spherical shape and can be calculated from the semi-principal axes of the ellipsoid used to approximate 1,3-DMU. The empirical hydrodynamic friction coefficient, C, linking macroscopic viscosity and molecular hydrodynamics, is usually found to be between its theoretical limits for stick, C_stick = 1, and slip, C_slip = 1 − f_⊥^−2/3, boundary conditions but for molecular reorientation is expected to be close to the slip limit. The present data yielded C = 0.25, which is significantly larger than C_slip = 0.129. This indicates strong interactions between solute and solvent and/or self-aggregation of 1,3-DMU in water.

More quantitative information was obtained from the relaxation strength of 1,3-DMU, S₁, which increased linearly when plotted against the molar concentration, c, of the solute (Fig. 3). This allowed for rearranging eqn (3) to yield

With the cavity factor of A₁ = 0.223, calculated from the semi-principal axes, μ_eff,1 = (11.3 ± 0.6) D was obtained for the effective dipole moment of 1,3-DMU. This value is significantly larger than permanent dipole moments, μ₁, determined for 1,3-DMU at infinite dilution in benzene (3.4 D), chloroform (4.4 D) or dichloroethane (3.8 D).⁷⁰ It also significantly exceeds the permanent moment of 3.8 D from quantum-chemical calculations and the corresponding apparent moment of 5.5 D in water (Table 2). This discrepancy indicates the pronounced parallel alignment of dipoles. Possible reasons for that could be the formation of hydrates where bound water molecules tightly follow solute dynamics, the formation of well-defined (i.e. long-lived) hydrates of 1,3-DMU molecules moving as a rigid entity or a more-or-less parallel but lose alignment of solute molecules, e.g. through dipole–dipole interactions or moderately strong hydrogen bonds between solute molecules (i.e. an effective Kirkwood factor g₁ > 1).

Table 2 Apparent dipole moments, μ_ap/D, of 1,3-DMU and its mono-, di- and tri-hydrates in vacuo and embedded in water estimated from quantum chemical calculations at the B3LYP/6-311G++(d,p) level of theory

Species	Species (Fig. 5)	μ
		Vacuum	Water
1,3-DMU	A	3.8	5.5
1,3-DMU·1H2O	B	4.9	7.2
	C	7.0	9.0
1,3-DMU·2H2O	D	6.2	10.6
	E	8.1	10.5
1,3-DMU·3H2O	F	8.9	13.9

---

An argument against the formation of stable solute aggregates is that over the entire concentration range the solute mode is well described by a Debye equation and the rotational correlation times calculated from τ₁ are linear in η over the entire concentration range (Fig. S2 of the ESI†). Thus, the existence of separate modes for monomeric and aggregated 1,3-DMU, a prerequisite for an aggregation equilibrium between two dipolar species, is unlikely. On the other hand, evaluation of the solvent relaxations indicated an equivalent of ∼1.6 water dipoles per solute molecule to be strongly bound in a hydrate. Additionally, ab initio quantum chemical calculations of 1,3-DMU·H₂O, 1,3-DMU·2H₂O and 1,3-DMU·3H₂O complexes at the B3LYP/6-311++G(d,p) level of theory, either in a vacuum or embedded in water, suggested that some of them have large dipole moments (Fig. 5, Table 2). As discussed in detail below, also the MD simulations support the formation of hydrates with coupled dynamics of 1,3-DMU and strongly interacting water molecules but there are also indications from MD for g₁ > 1 through moderately strong 1,3-DMU–1,3-DMU hydrogen bonds. As discussed in detail below, combining the results of the present quantum chemical calculations with the analysis of contact and H-bond lifetimes times from MD allows us to suggest that 2H₂O molecules H-bonded to the solute's carbonyl oxygen are most likely responsible for the strongly bound water fraction observed in DRS. Indeed, the dipole moment of the corresponding hydrate (structure D of Fig. 5) is in good agreement with the experimental value of μ_eff,1 = (11.3 ± 0.6) D.

Fig. 5 Structures of 1,3-DMU and its mono-, di- and tri-hydrates from DFT calculations along with their dipole moment vectors.

Molecular dynamics simulations

1,3-DMU hydration and aggregation. Much knowledge about solute hydration can be obtained from studying radial distribution functions (RDFs). These allow, among others, us to calculate the number of water or other solute molecules in direct contact with a solute molecule, i.e. determination of the first-shell coordination numbers CN_MD(H₂O) and CN_MD(DMU). We calculated these numbers from proximal radial distribution functions (pRDFs)^71,72 for the solute–solvent and solute–solute distributions in the 1,3-DMU–water mixtures.

The pRDFs for the 1,3-DMU–O_W and 1,3-DMU–C pairs are plotted in Fig. 6 (see the inset for the atom labeling on 1,3-DMU). Note that these pRDFs are not normalized by 4πr². As the solute concentration does not greatly affect the positions of the maxima and the minima of these functions, but only their height or depth, only the pRDFs obtained at c = 3.0 M are shown (see Table S5 of the ESI† for more information). As can be seen, there are two peaks at ∼2.0 and ∼2.7 Å on the 1,3-DMU–O_W pRDF. The first small peak is due to the interactions of the H_N and H_C atoms of 1,3-DMU with the oxygen atom of water; accordingly, the corresponding minimum at 2.2 Å does not correspond to the solute's first hydration shell (FHS). The second maximum indicates the interactions of water with the heavy solute atoms (N, O, C, C_H). The corresponding minimum observed at 3.8 Å defines the complete first hydration shell of the solute comprising 25.9H₂O molecules at infinite dilution. With increasing concentration CN_MD(H₂O) continuously drops, reaching a value of 12.5 at c = 4.9 M (Fig. 4).

Fig. 6 Proximal radial distribution functions (pRDFs) of water oxygen atoms (O_W) and 1,3-DMU carbons (C) relative to the atoms of a reference solute obtained at 3.0 M. Arrows address the curves to the appropriate axes. The inset displays the atom labeling used for 1,3-DMU in this work.

Comparison of CN_MD(H₂O) with the total effective hydration number from DRS (Fig. 4) yields Z_t > CN_MD(H₂O) for c < 1.5 M. This may indicate that at low solute concentrations also H₂O molecules beyond the first hydration shell of 1,3-DMU are retarded in their dynamics. Currently, this cannot be confirmed by other methods. However, the present MD simulations clearly show (see below) that only a small fraction of the solvent molecules residing in the first hydration shell, equivalent to Z_ib ≈ 1.6 frozen H₂O dipoles, is impeded so strongly in its rotational dynamics that these molecules are neither detectable as bulk-like nor as slow water by DRS.

Note that CN_MD(H₂O) > Z_t at c > 2.0 M (Fig. 4) and eventually exceeding the solvent-to-solute ratio is a consequence of the strong solvation shell overlap, which is not taken into account in the calculation of CN_MD(H₂O). Instead, all H₂O molecules surrounding a given 1,3-DMU molecule are counted irrespective if they are just in contact with this reference or are shared with neighbouring solute molecules.

Direct comparison of the present values from MD and DRS for 1,3-DMU with effective hydration numbers from compressibility studies, namely 4.7⁷³ and 4.3⁷⁴ at c → 0, is not straightforward but obviously also this technique monitors only a part of the first hydration shell. On the other hand the NMR study of Costantino et al.⁵⁰ yielded a constant number, 19.79 ± 0.27, over the concentration range 0.061 ≤ c/M ≤ 1.76. According to the authors this number encompasses water molecules solvating both hydrophilic and hydrophobic moieties of 1,3-DMU and thus should correspond to Z_t= Z_ib + Z_s from DRS.⁷⁵ However, since the NMR values are comparable to the first-shell coordination numbers, CN_MD(H₂O), but smaller than Z_t for c < 1.7 M this indicates a “view range” of NMR restricted to the first hydration shell.

The 1,3-DMU–C pRDF describing solute–solute contacts exhibits a well-defined structure (Fig. 6). Two peaks at 3.5 and 5.1 Å are observed for the 3.0 M solution. Here each 1,3-DMU molecule is in contact with, on average, 1.7 other 1,3-DMU molecules within the radius of 4.1 Å defined by the first minimum. This number increases with increasing concentration from CN_MD(DMU) = 0.3 at 0.48 M to CN_MD(DMU) = 2.5 at 4.9 M. Inclusion of the second peak, i.e. integrating the pRDF up to 5.7 Å, gives 3.8 neighboring 1,3-DMU molecules at 0.48 M while at the highest studied concentration each 1,3-DMU is surrounded by 5.9 solute molecules within this distance. Similar to the 1,3-DMU–O_W pRDF, the pronounced shoulder on the 1,3-DMU–C pRDF's at ∼3.0 Å is caused by H_N and H_C atoms interacting with the C atom of the reference molecule.

Further details on solute–solvent and solute–solute interactions can be extracted from the RDFs for individual atom pairs on solvent and solute. As the latter has two possible centers able to form hydrogen bonds, namely the carbonyl oxygen atom as an acceptor and the amino hydrogen atoms as donors, the most interesting RDFs are O–H_W, H_N–O_W, and H_N–O correlations. For the 3.0 M solution these are shown in Fig. 7.

Fig. 7 Radial distribution functions (RDFs) between H₂O and 1,3-DMU atoms in the 3.0 M solution.

In the O–H_W RDF, the first peak corresponds to the direct coordination of the carbonyl oxygen to the hydrogen atoms of water. The second peak represents water molecules with an H_W atom not directly bound to the O atom of 1,3-DMU. The third broad peak at about 5.2 Å corresponds to H₂O molecules coordinated to the remaining part of the solute and to H₂O molecules belonging to the second hydration shell.

With regard to the amino groups, the first peak in the H_N–O_W RDF is located at 2.0 Å. At about 5.2 Å, a second – rather broad – peak corresponds to H₂O molecules located either in the second hydration shell or close to the carbonyl O. Interestingly, for the H_N–O_W and H_N–O RDFs the positions of their first maxima practically coincide. This suggests that – at least at high concentrations – not only H-bonded hydrates are formed but also head-to-tail aggregates of 1,3-DMU molecules via amino groups interacting with carbonyl oxygen.

A detailed picture of the arrangement of the water molecules hydrogen bonded to 1,3-DMU is given by the spatial distribution function of the three closest water molecules, shown in Fig. 8. There are indeed two water molecules interacting with the carbonyl oxygen atom acting as the H-bond donor. Additionally, a further H₂O molecule interacts simultaneously with both NH acceptor groups. The two water molecules bound to the CO group show up as bent toroidal clouds, the water molecule that is hydrogen bonded to the NH moieties forms a nearly ellipsoidal cloud located directly above the mid-point between these groups. The analysis of combined distribution functions linking the distance between the C atom of 1,3-DMU and water oxygen, O_W, and relative orientations of water molecules hydrating the solute provided in the ESI† reveals that the orientations found for the closest three water molecules are consistent with the picture that they form hydrogen bonds to 1,3-DMU and to a certain extent align their dipole vectors with that of the solute. This picture agrees qualitatively with the minimum-energy structures of 1,3-DMU–H₂O complexes shown in Fig. 5 and explains the experimental result of Z_ib ≈ 1.6 “frozen” water dipoles per solute molecule and the large effective dipole moment of μ_eff,1 = (11.3 ± 0.6) D of 1,3-DMU.

Fig. 8 Spatial distribution functions of the atoms of the nearest three H₂O molecules around 1,3-DMU in the 0.48 M solution.

Fig. 9 shows the spatial distribution functions of the atoms of neighbouring 1,3-DMU molecules around the solute. Clearly, the bent oblate ellipsoid zones formed by the O, H_N, C and N atoms above and below the CNN plane of the reference 1,3-DMU molecule indicate the tendency of solute molecules to form stacked structures. Simultaneously, the spherical distribution of O atoms just on top of the amino hydrogens of the reference, as well as the regions of H atoms in the vicinity of the carbonyl oxygen, suggests the formation of head-to-tail associates stabilized by H-bonding. At c = 3.0 M the solute dipole vectors are aligned roughly parallel with respect to each other in ∼36% of these structures, whereas in ∼23% they are anti-parallel. A more detailed analysis of the arrangement of nearest 1,3-DMU molecules around the reference solute, performed by means of a 2D analysis combining radial and angular distribution functions, is given in the ESI.†

Fig. 9 Spatial distribution functions of selected atoms of 1,3-DMU molecules closest to the reference solute molecule obtained for the 3.0 M solution.

1,3-DMU reorientation. Reorientational motions of molecules in liquids can be analyzed through time correlation functions

C_u,l(t) = 〈P_l[u(t) × u(0)]〉, (10)

of appropriate vectors, u, where P_l(x) is a Legendre polynomial of degree l.⁷⁶ For dielectric relaxation, monitoring the reorientation of the molecular dipole vector, l = 1. Experiments sensitive to tensorial properties, like NMR or time-resolved infrared spectroscopy, probe l = 2 whereas the reorientation corresponding to l = 3 is currently available only through MD simulations.⁷⁷ The behavior of C_u,1(t), C_u,2(t) and C_u,3(t) can be characterized by rotational correlation times τ_u,1, τ_u,2, and τ_u,3, respectively, which are defined as the time integral from zero to infinity of C_u,l(t).

The correlation times for the reorientation of the three vectors depicted in the inset of Fig. 10 were derived from the MD simulations. The first vector, x, points along the C=O bond and coincides with the dipole vector of 1,3-DMU. The second vector, y, which points along the N–N axis, lies in the molecular plane and is perpendicular to the first. The third vector, z, is perpendicular to both of them and thus normal to the molecular plane. To calculate their correlation times the corresponding C_u,l(t) functions were fitted to exponential decays and integrated analytically. The obtained values for τ_u,l are listed in Table 3.

Fig. 10 Rotational correlation times from DRS, τ_rot,1, as a function of solute concentration and corresponding data, τ_x,1′, τ_y,1′, and τ_z,1′, for the reorientation of the x, y, and z-vectors describing 1,3-DMU orientation as obtained by MD simulations.

Table 3 First, τ₁, second, τ₂, and third-order, τ₃, rotational correlation times of the vectors x, y and z describing the orientation of 1,3-DMU obtained from MD simulations

c/M	x			y			z
	τ 1	τ 2	τ 3	τ 1	τ 2	τ 3	τ 1	τ 2	τ 3
0.4755	18.3	6.5	1.7	26.6	8.7	2.3	14.1	4.4	0.9
0.9327	21.0	8.1	1.9	31.5	9.9	2.6	16.1	5.1	1.0
1.7373	30.8	11.0	2.5	33.8	11.1	2.9	21.7	6.8	1.2
3.0463	46.3	16.3	3.6	43.5	15.3	3.8	30.3	9.9	1.6
4.0780	64.5	21.7	4.8	49.1	17.4	4.5	38.8	12.3	1.9
4.9172	94.0	30.7	6.4	59.9	22.5	5.4	47.1	15.4	2.4

---

In the case of molecular reorientation through rotational diffusion the ratio τ_u,1/τ_u,2 should be equal to 3. Inspection of the τ_u,1 and τ_u,2 data from MD summarized in Table 3 reveals 2.7 ≤ τ_u,1/τ_u,2 ≤ 3.1, depending on the chosen vector and solute concentration. This would suggest 1,3-DMU reorientation through rotational diffusion. However, as shown by Laage et al.,⁷⁷ comparison of first- and second-order reorientation times may not be sufficient to identify the relaxation mechanism and additionally the ratio τ_u,1/τ_u,3 has to be considered. For rotational diffusion theory predicts τ_u,1/τ_u,3 = 6. In the present case the simulations yielded for the 3.0 M solution the ratios of 12.9, 11.4 and 18.8 for the vectors x, y and z, suggesting that the reorientation of 1,3-DMU molecules cannot be classified as rotational diffusion.

Due to missing far-infrared/terahertz data and the decomposition of the spectra into individual modes according to eqn (2) DRS does not yield the time integral of the collective dipole correlation function. Instead, the dielectric relaxation time of 1,3-DMU, τ₁ (Table S2 of the ESI†), is given by the slope of this function at long times. The corresponding rotational correlation time, τ_rot,1, of a 1,3-DMU molecule, obtained from τ₁ with the Powles–Glarum equation,^67,68 is compared in Fig. 10 with the long-time slope, τ_u,1′, of the rotational correlation function, C_u,1(t), from the MD simulations (see Table S8 of the ESI† for details).

The very good general agreement between the experimental data and the rotational correlation times from MD immediately confirms our assignment of the lowest frequency mode (amplitude S₁, relaxation time τ₁) resolved in the dielectric spectra to the reorientation of 1,3-DMU dipoles. At the same token, this good agreement between experiment and simulation validates the force field used in this MD study. The combination of the SPC/E model for water and a customized KBFF model for 1,3-DMU is clearly able to reproduce the dynamic properties of 1,3-DMU–water mixtures.

Up to c ≈ 2.0 M the values for τ_x,1′, τ_y,1′, and τ_z,1′ remain rather similar. However, at higher concentrations the reorientation of the vector parallel to the C=O bond is considerably slowed down, i.e. τ_x,1′ increases significantly. This suggests that, similar to urea,^29,78–80 the reorientation of 1,3-DMU dipoles is becoming more and more anisotropic and possibly indicates increasing 1,3-DMU aggregation in line with the structural information from MD (Fig. 9) and the inference of Rodríguez-Ropero et al.²³

Contact times and H-bond dynamics. The dielectric relaxation times discussed previously revealed a significant slowdown of water dynamics induced by the presence of 1,3-DMU. The fractions of slow (s) and bound (ib), i.e. of moderately and strongly retarded water observed in DRS (Fig. 4), are almost certainly caused by specific interactions of H₂O molecules with different moieties of the 1,3-DMU molecule. Estimation of the rotational correlation time of the solvent in the vicinity of specific entities of the solute is, however, hampered by computational difficulties associated with the rather fast solvent exchange (see below) that causes poor statistics at long times. On the other hand, due to the collective wait, jump and settle mechanism of water reorientation, which also applies to the H₂O molecules in the hydration shell, albeit with waiting times different from the bulk, solvent exchange will also influence the experimental τ₂ values. Thus, in order to reveal the connection between water dynamics and interaction sites we analyzed the contact times, τ_C, of H₂O molecules with hydrophobic and hydrophilic moieties of the solute. To estimate the dynamics of the H-bonds between appropriate particles/sites we also calculated corresponding H-bond lifetimes, τ_H. As urea derivatives are prone to self-aggregation,^70,81–84 the contact times and H-bond lifetimes between neighboring 1,3-DMU molecules were also calculated. To do so, we define an aggregate existence correlation function C(t) as

C(t) = 〈a_ij(0)a_ij(t)〉/〈a_ij(0)²〉, (11)

where a_ij is the value of a binary parameter for a pair of sites i and j. The variable a_ij equals one if the declared conditions for the observed sites are satisfied at the given moment, t, and zero otherwise. Such a definition corresponds to the so-called intermittent correlation function as it allows two molecules to dissociate and re-associate during the period between time zero and t. After correcting for finite size effects the corresponding lifetimes were obtained by the integrationIn the case of contact times, only a radial cut-off was used as a geometric criterion. The respective cut-off distances were taken as the first minima of the individual RDFs for each particular site – O_W (H₂O) or site – C (1,3-DMU) interaction. For H-bonding, a combination of radial and angular cut-offs was applied. The corresponding RDFs and combined radial/angular distribution functions (CDFs) for these interactions are shown in Fig. S8 and S9 of the ESI,† and the values of the corresponding radial, r^c, and angular, θ^c, cut-offs as well as the contact and H-bond lifetimes are listed in Table 4.

Table 4 Calculated intermittent contact times, τ_C, and H-bond lifetimes, τ_H, for O_W–O_W, O–O_W, H_N–O_W, C_H–O_W, and C–C, H_N–O interactions as well as the corresponding radial, r^c, and angular, θ^c, cut-offs

c/M	τ C/ps					τ H/ps
	OW–OW	O–OW	HN–OW	CH–OW	C–C	OW–OW	O–OW	HN–OW	HN–O
	r ^c = 3.3 Å	r ^c = 3.3 Å	r ^c = 2.7 Å	r ^c = 5.1 Å	r ^c = 5.0 Å	r ^c = 3.3 Å, θ^c = 30°	r ^c = 3.3 Å, θ^c = 30°	r ^c = 3.3 Å, θ^c = 40°	r ^c = 3.4 Å, θ^c = 40°
0.0000	7.9					7.1
0.4755	8.6	11.2	5.1	20.6	55.1	7.8	10.7	4.1	9.2
0.9327	9.3	12.6	5.8	22.1	51.0	8.5	12.1	4.7	12.5
1.7373	10.7	16.0	7.4	24.9	49.1	9.8	15.4	6.0	21.9
3.0463	13.7	22.1	10.6	30.4	57.6	12.5	21.3	8.9	34.0
4.0780	17.0	28.1	14.0	35.8	64.7	15.5	27.1	11.9	45.7
4.9172	19.8	36.5	19.1	42.2	73.8	18.2	35.4	16.5	60.3

---

As seen from Table 4, the contact times of H₂O molecules in the vicinity of the various 1,3-DMU moieties, respectively, and H₂O itself, decrease in the sequence C_H > O > H₂O > H_N for all 1,3-DMU–H₂O mixtures studied. Compared to water surrounded by water, the exchange of H₂O near the O atom of the solute is retarded by a factor of 1.6, whereas H₂O molecules near the amino groups exchange 1.3 times faster for the 3.0 M solution. This result is in line with an earlier MD investigation of urea–water mixtures by Kokubo et al.⁸⁵ using the KBFF model for urea and TIP3P water. The authors showed that the H_N atoms of urea are less effective hydrogen-bond donors than water itself due to the pronounced nucleophilicity of the carbonyl group. In the case of 1,3-DMU, methylation of the nitrogen atoms leads to an additional increase in the electron density of the oxygen lone pair orbitals and a decrease in the partial electron density of the nitrogen lone pair orbitals.⁸⁶ This manifests in a less negative partial charge for nitrogen in 1,3-DMU (q_N = −0.553) compared to urea (q_N = −0.693⁴⁵). Thus, methyl substitution reduces solvent interaction with the nitrogen atom and promotes it at the oxygen atom.

Comparison of τ_C and τ_H revealed that, similar to water surrounded by water, the contact times for H₂O around the carbonyl and amino groups of the solute are mainly determined by H-bonding. For example, at 3.0 M the values of τ_C = 22.1 ps and τ_H = 21.3 ps were found for O–H₂O interactions (Table 4). The corresponding data for H₂O–H₂O are τ_C = 13.7 ps and τ_H = 12.5 ps. Thus, O–H₂O dynamics is slowed down by a factor of ∼1.6 compared that of H₂O surrounded by H₂O. Surprisingly, however, the dynamics of the H_N–H₂O interactions is notably faster with τ_C = 10.6 ps and τ_H = 8.9 ps. On the other hand, as predicted by the excluded volume mechanism in the study of Laage et al.,^60,87 the contact time of H₂O molecules with the methyl groups of 1,3-DMU, τ_C = 30.4 ps, is significantly (2.2 fold) increased.

When comparing these contact times with experimental dielectric relaxation times it has to be considered that the simulation data do not directly reflect the reorientation of the water-dipole vector. On the other hand, water reorientation is not diffusive but, due to H-bonding, occurs through sub-picosecond large-angle jumps separated by waiting times of several picoseconds.^60,87 Whilst the sub-picosecond jumps can be associated with the present fast water mode (S₄, τ₄),⁵⁸ the environment-dependent waiting time between these jumps determines the dielectric relaxation times τ₂ and τ₃. Although numerical equivalence cannot be expected, it is reasonable to assume that τ₃ should be comparable to τ_C of H₂O surrounded by H₂O. Similarly, τ₂, which is characteristic of hydration water impeded in its dynamics, should have an equivalent in the contact times of hydration water with appropriate solute moieties.

The experimental retardation factor (τ₂/τ₃ = 2.0 at 3.0 M) is indeed comparable to the relative increase of H₂O contact times around the carbonyl (τ_C(O–O_W)/τ_C(O_W–O_W) = 1.6) and the methyl groups (τ_C(C_H–O_W)/τ_C(O_W–O_W) = 2.2) of 1,3-DMU. At least for c ≥ 1.5 M where Z_s ≤ CN_MD and solute hydration shells start to overlap (Fig. 4), we may therefore assign the DRS-detected slow water fraction to H₂O molecules in contact with these two groups of 1,3-DMU.

Note that from τ_C(H_N–O_W)/τ_C(O_W–O_W) < 1 a decrease of the dielectric relaxation time τ₃ could be expected or even the emergence of a separate mode for water hydrating the amino groups. However, according to the present simulations only a single H₂O molecule (the 3rd-next neighbor) shared by both amino groups is involved here (Fig. 8; also see the discussion in the ESI†). The contribution of this molecule is certainly too small to give rise to a detectable DRS mode. However, due to a likely band overlap it might be a possible reason for the rather strong increase of the fast-water relaxation time, τ₄, with concentration (Table S2, ESI†).

The present MD simulations revealed that the dipole moment vectors of the three H₂O molecules closest to the solute molecule, with neighbors 1 & 2 hydrogen bonding to the carbonyl group and the third shared by the amino groups essentially, wobble in a 30° cone around the dipole direction (x-vector) of 1,3-DMU (Fig. 10; also see Fig. 8 and the discussion in the ESI;† note that this structure is rather similar to hydrate cluster F in Fig. 5). Thus, their components normal to x will either contribute to Z_s (neighbors 1 & 2) or to the bulk-water relaxation (neighbor 3) whereas the “frozen” parallel components give rise to Z_ib ≈ 1.6 and the large effective dipole moment, μ_eff,1 = (11.3 ± 0.6) D. Despite its short contact time also neighbor 3 may be involved here as the H_N–O_W RDF exhibits a well-defined maximum (Fig. S8, ESI†) clearly associated with hydrogen bonding (Fig. S9, ESI†).

It is noteworthy that the contact and H-bond lifetimes of direct 1,3-DMU–1,3-DMU interactions summarized in Table 4 are significantly longer than those for H₂O molecules both near H₂O and the solute but comparable to the rotational correlation time of the solute found with DRS (Fig. 10). According to the MD simulations solute–solute aggregates (Fig. 9 and Fig. S6, ESI†) live up to 57.6 ps and can remain H-bonded for about 34 ps at 3.0 M (Table 4), whereas the correlation time for the reorientation of the CNN plane of a 1,3-DMU molecule is τ_y,1′ = 59.0 ps (Table S8, ESI†). These values are similar to the experimental dielectric relaxation time, τ₁ = 67 ps (Table S2, ESI†), of the solute at this concentration and suggest that the formation of such aggregates is a likely reason for the considerable increase of τ₁ (and possibly also of viscosity) with solute concentration. The analysis of the number, n_H, of H-bonds formed between solute molecules (Fig. S11, ESI†) showed indeed that the fraction of fully hydrated 1,3-DMU molecules (n_H = 0) decreases from (98.5 ± 1.5)% at 0.48 M to (67.6 ± 5.2)% at 4.9 M. Therefore, in the most concentrated solution studied ∼30% of the solute molecules form head-to-tail associates via one (8.0 ± 0.6)% or even two (22.3 ± 4.3)% hydrogen bonds. Simultaneously, the tendency of anti-parallel stacking increases, see the detailed discussion in the ESI.† As a consequence, solute hydration shells increasingly overlap, leading to the pronounced decrease of the effective hydration number, Z_s, associated with “slow” water from ∼39 at c → 0 to 3.3 at c = 4.9 M. Since the number of H₂O molecules hydrogen bonded to the carbonyl groups of 1,3-DMU decreases only slightly (Fig. S11, ESI†) the number of “frozen” water dipoles remains practically constant at Z_ib ≈ 1.6. The opposing effects on the effective dipole moment of 1,3-DMU arising from the simultaneous formation of parallel head-to-tail aggregates and anti-parallel stacks apparently compensate, so that μ_eff,1 = (11.3 ± 0.6) D does not change with solute concentration.

Conclusions

In this contribution we studied 1,3-DMU–water mixtures combining dielectric relaxation spectroscopy (DRS) and molecular dynamics (MD) simulations. We found that the static permittivity increases considerably in solutions of 1,3-DMU relative to pure water due to the large effective dipole moment of the solute. The latter arises from the pronounced parallel alignment of the dipole moments of 1,3-DMU and its three closest H₂O molecules, of which two form long-lived hydrogen bonds to the carbonyl oxygen atom of the solute, whereas the third forms surprisingly short-lived H-bonds to both amino groups. This configuration leads to Z_ib ≈ 1.6 “frozen” water dipoles independent of concentration. With increasing concentration solute aggregation via anti-parallel stacking and head-to-tail H-bonding becomes relevant. The associated overlap of hydration shells is reflected by a strong decrease of the number of “slow” (retarded) water molecules, Z_s, hydrating 1,3-DMU. Interestingly, due to compensation effects this aggregate formation does neither affect Z_ib nor μ_eff,1.

Acknowledgements

V. A. acknowledges support from the Russian Government Program for Competitive Growth of Kazan Federal University and thanks the Deutscher Akademischer Austauschdienst (DAAD) for a grant enabling her experiments in Regensburg.

Notes and references

1. E. Freire, A. Schön, B. M. Hutchins and R. K. Brown, Drug Discovery Today, 2013, 18, 1007–1013.
2. R. Kumaran and P. Ramamurthy, J. Fluoresc., 2011, 21, 1499–1508.
3. D. R. Canchi and A. E. García, Annu. Rev. Phys. Chem., 2013, 64, 273–293.
4. D. Horinek and R. R. Netz, J. Phys. Chem. A, 2011, 115, 6125–6136.
5. B. Moeser and D. Horinek, J. Phys. Chem. B, 2014, 118, 107–114.
6. W. Zheng, A. Borgia, K. Buholzer, A. Grishaev, B. Schuler and R. B. Best, J. Am. Chem. Soc., 2016, 138, 11702–11713.
7. T. T. Herskovits, H. Jaillet and B. Gadegbeku, J. Biol. Chem., 1970, 245, 4544–4550.
8. A. Shimizu, K. Fumino, K. Yukiyasu and Y. Taniguchi, J. Mol. Liq., 2000, 85, 269–278.
9. K. Hamaguchi and H. Sakai, J. Biochem., 1965, 57, 721–732.
10. J. Quijada, G. López, R. Versace, L. Ramírez and M. L. Tasayco, Biophys. Chem., 2007, 129, 242–250.
11. V. M. Tischenko, Biophysics, 2011, 56, 602–605.
12. R. Zangi, R. Zhou and B. J. Berne, J. Am. Chem. Soc., 2009, 131, 1535–1541.
13. D. R. Canchi, D. Paschek and A. E. García, J. Am. Chem. Soc., 2010, 132, 2338–2344.
14. S. Funkner, M. Havenith and G. Schwaab, J. Phys. Chem. B, 2012, 116, 13374–13380.
15. K.-J. Tielrooij, J. Hunger, R. Buchner, M. Bonn and H. J. Bakker, J. Am. Chem. Soc., 2010, 132, 15671–15678.
16. Y. L. A. Rezus and H. J. Bakker, Phys. Rev. Lett., 2007, 99, 148301.
17. Y. L. A. Rezus and H. J. Bakker, J. Phys. Chem. A, 2008, 112, 2355–2361.
18. J. Hunger, K.-J. Tielrooij, R. Buchner, M. Bonn and H. J. Bakker, J. Phys. Chem. B, 2012, 116, 4783–4795.
19. N. Poklar, G. Vesnaver and S. Lapanje, Biophys. Chem., 1993, 47, 143–151.
20. N. Poklar, G. Vesnaver and S. Lapanje, J. Protein Chem., 1995, 14, 709–719.
21. N. Poklar, G. Vesnaver and S. Lapanje, Biophys. Chem., 1996, 57, 279–289.
22. L. B. Sagle, Y. Zhang, V. A. Litosh, X. Chen, Y. Cho and P. S. Cremer, J. Am. Chem. Soc., 2009, 131, 9304–9310.
23. F. Rodríguez-Ropero and N. F. A. van der Vegt, J. Phys. Chem. B, 2014, 118, 7327–7334.
24. S. Subramanian, T. S. Sarma, D. Balasubramanian and J. C. Ahluwalia, J. Phys. Chem., 1971, 75, 815–820.
25. U. Kaatze, H. Gerke and R. Pottel, J. Phys. Chem., 1986, 90, 5464–5469.
26. P. R. Philip, G. Perron and J. E. Desnoyers, Can. J. Chem., 1974, 52, 1709–1713.
27. E. V. Ivanov and V. K. Abrosimov, J. Chem. Eng. Data, 2013, 58, 1103–1111.
28. R. Buchner and G. Hefter, Phys. Chem. Chem. Phys., 2009, 11, 8984–8999.
29. V. Agieienko and R. Buchner, Phys. Chem. Chem. Phys., 2016, 18, 2597–2607.
30. T. Shikata and S. Itatani, J. Solution Chem., 2002, 31, 823–844.
31. A. Eiberweiser, A. Nazet, S. E. Kruchinin, M. V. Fedotova and R. Buchner, J. Phys. Chem. B, 2015, 119, 15203–15211.
32. A. Schönhals and F. Kremer, in Broadband Dielectric Spectroscopy, ed. F. Kremer and A. Schönhals, Springer Berlin Heidelberg, 2003, ch. 1, pp. 1–33.
33. T. Sonnleitner, D. A. Turton, S. Waselikowski, J. Hunger, A. Stoppa, M. Walther, K. Wynne and R. Buchner, J. Mol. Liq., 2014, 192, 19–25.
34. J. Barthel, R. Buchner, P. N. Eberspächer, M. Münsterer, J. Stauber and B. Wurm, J. Mol. Liq., 1998, 78, 83–109.
35. R. Buchner, T. Chen and G. Hefter, J. Phys. Chem. B, 2004, 108, 2365–2375.
36. A. Stoppa, A. Nazet, R. Buchner, A. Thoman and M. Walther, J. Mol. Liq., 2015, 212, 963–968.
37. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford CT, 2004.
38. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
39. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
40. J. B. Foresman, T. A. Keith, K. B. Wiberg, J. Snoonian and M. J. Frisch, J. Phys. Chem., 1996, 100, 16098–16104.
41. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129–145.
42. B. H. Besler, K. M. Merz and P. A. Kollman, J. Comput. Chem., 1990, 11, 431–439.
43. C. Oostenbrink, A. Villa, A. E. Mark and W. F. Van Gunsteren, J. Comput. Chem., 2004, 25, 1656–1676.
44. S. Weerasinghe and P. E. Smith, J. Chem. Phys., 2003, 118, 5901–5910.
45. S. Weerasinghe and P. E. Smith, J. Phys. Chem. B, 2003, 107, 3891–3898.
46. K. J. Miller and J. Savchik, J. Am. Chem. Soc., 1979, 101, 7206–7213.
47. J. Hinze and H. H. Jaffe, J. Am. Chem. Soc., 1962, 84, 540–546.
48. M. G. Bravo-Sánchez, A. Estrada-Baltazar, G. A. Iglesias-Silva and M. Ramos-Estrada, J. Chem. Eng. Data, 2009, 55, 989–991.
49. R. Sartorio, P. Padulano, L. Costantino and V. Vitagliano, J. Solution Chem., 1981, 10, 111–119.
50. L. Costantino, G. D'Errico, O. Ortona and V. Vitagliano, J. Mol. Liq., 2000, 84, 179–191.
51. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
52. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
53. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
54. S. Nosé and M. L. Klein, Mol. Phys., 1983, 50, 1055–1076.
55. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
56. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
57. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.
58. T. Fukasawa, T. Sato, J. Watanabe, Y. Hama, W. Kunz and R. Buchner, Phys. Rev. Lett., 2005, 95, 197802.
59. A. Eiberweiser, A. Nazet, G. Hefter and R. Buchner, J. Phys. Chem. B, 2015, 119, 5270–5281.
60. D. Laage, G. Stirnemann, F. Sterpone, R. Rey and J. T. Hynes, Annu. Rev. Phys. Chem., 2011, 62, 395–416.
61. R. Buchner, C. Holzl, J. Stauber and J. Barthel, Phys. Chem. Chem. Phys., 2002, 4, 2169–2179.
62. H. M. A. Rahman, G. Hefter and R. Buchner, J. Phys. Chem. B, 2012, 116, 314–323.
63. R. Buchner, J. Barthel and J. Stauber, Chem. Phys. Lett., 1999, 306, 57–63.
64. J. Barthel, H. Hetzenauer and R. Buchner, Ber. Bunseng. Phys. Chem., 1992, 96, 1424–1432.
65. H. M. A. Rahman, G. Hefter and R. Buchner, J. Phys. Chem. B, 2013, 117, 2142–2152.
66. Note that the present CN numbers are calculated without taking into accout possible solvation shell overlap and associated sharing of H₂O molecules.
67. J. G. Powles, J. Chem. Phys., 1953, 21, 633–637.
68. S. H. Glarum, J. Chem. Phys., 1960, 33, 639–643.
69. J. L. Dote and D. Kivelson, J. Phys. Chem., 1983, 87, 3889–3893.
70. M. Obrzud, M. Rospenk and A. Koll, J. Phys. Chem. B, 2010, 114, 15905–15912.
71. P. K. Mehrotra and D. L. Beveridge, J. Am. Chem. Soc., 1980, 102, 4287–4294.
72. B. Lin and B. M. Pettitt, J. Chem. Phys., 2011, 134, 106101.
73. A. Burakowski and J. Gliński, Chem. Phys. Lett., 2015, 641, 40–43.
74. J. Krakowiak, J. Wawer and A. Panuszko, J. Chem. Thermodyn., 2013, 58, 211–220.
75. Costantino et al.⁵⁰ assume viscosity to be independent of the solute concentration. Correction of their diffusion data with the present η values (Table S1, ESI†) yields an even higher NMR hydration number of 22 ± 1.
76. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Elsevier Science, 2006.
77. D. Laage and J. T. Hynes, Science, 2006, 311, 832–835.
78. A. Idrissi, F. Sokolić and A. Perera, J. Chem. Phys., 2000, 112, 9479–9488.
79. A. Idrissi, P. Bartolini, M. Ricci and R. Righini, Phys. Chem. Chem. Phys., 2003, 5, 4666–4671.
80. P. O. Åstrand, A. Wallqvist, G. Karlström and P. Linse, J. Chem. Phys., 1991, 95, 8419–8429.
81. R. Adams, H. H. M. Balyuzi and R. E. Burge, J. Appl. Crystallogr., 1977, 10, 256–261.
82. J. Jadżyn, L. Bouteiller, J.-L. Déjardinc and G. Czechowskia, Acta Phys. Pol., A, 2006, 110, 495–505.
83. M. Obrzud, M. Rospenk and A. Koll, J. Mol. Struct., 2012, 1018, 54–63.
84. J. Świergiel, L. Bouteiller and J. Jadżyn, J. Phys. Chem. B, 2015, 119, 12947–12953.
85. H. Kokubo and B. M. Pettitt, J. Phys. Chem. B, 2007, 111, 5233–5242.
86. E. I. Harnagea and P. W. Jagodzinski, Vib. Spectrosc., 1996, 10, 169–175.
87. D. Laage, G. Stirnemann and J. T. Hynes, J. Phys. Chem. B, 2009, 113, 2428–2435.

---

Footnote

† Electronic supplementary information (ESI) available: Tables with data for density, viscosity and with relaxation parameters describing the dielectric spectra; tables with bonded and non-bonded parameters for the solute force field; auxiliary tables and figures. See DOI: 10.1039/c6cp07407c

---