Hydrogen-bond dynamics of water in presence of an amphiphile, tetramethylurea: signature of confinement-induced effects

Various experimental and simulation studies have suggested that the presence of amphiphilic molecules in aqueous solutions substantially perturbs the tetrahedral hydrogen-bond (H-bond) network of neat liquid water. Such structural perturbation is expected to impact H-bond lifetime of liquid water. Tetramethylurea (TMU) is an example of an amphiphile because it possesses both hydrophobic and hydrophilic moieties. Molecular dynamics simulations of (water+TMU) binary mixtures at various compositions have been performed in order to investigate the microscopic mechanism through which the amphiphiles influence the H-bond dynamics of liquid water at room temperature. Present simulations indicate lengthening of both water–water H-bond lifetime and H-bond structural relaxation time upon addition of TMU in aqueous solution. At the highest TMU mole fraction studied, H-bond lifetime and structural relaxation time are, respectively, ∼4 and ∼8 times longer than those in neat water. This is comparable with the slowing down of H-bond dynamics for water molecules confined in cyclodextrin cavities. Simulated relaxation profiles are multi-exponential in character at all mixture compositions, and simulated radial distribution functions suggest enhanced water–water and water–TMU interactions upon addition of TMU. No evidence for complete encapsulation of TMU by water H-bond network has been found.

Behaviour of water molecules in the vicinity of hydrophobic solute molecules is an issue of intense research. It has been found from both experimental and theoretical studies that reorientation dynamics of water molecules in the hydration shell of hydrophobes exhibits substantial retardation. [36][37][38][39][40][41][42][43] Dynamic anisotropy measurements of water around a hydrophobic molecule [32] reveal that water molecules become immobilised in the hydration shell of hydrophobic solute, supporting the controversial iceberg model. [44] Interestingly, a combined computer simulation and theoretical study [38] has suggested no immobilisation of water molecules in the hydrophobic hydration layer and thus debated the validity of the iceberg model. Soluteexcluded volume effects at the transition state are found to be responsible for the slowing down of water reorientation. Concentrated aqueous solutions are found to behave like confining environments due to the steric effects rendered by the solute molecules. [45,46] Hydrophobic molecules tend to aggregate in solution. At higher solute concentrations, these aggregates distribute throughout the solution. Consequently, water molecules are forced to reside by the side of the surfaces created by these aggregated moieties. This results in significant retardation of reorientational and translational dynamics of interfacial water molecules residing at the surface of the solute aggregates. [47] By this way, water shows similar behaviour in concentrated aqueous solutions as found in explicit confinements.
The extensive hydrogen-bond (H-bond) network of water plays a key role in determining the structure and dynamics of aqueous systems. Therefore, experiments, [48 -52] theory and simulations [53 -60] have been employed to gain insight into the mechanism of both H-bond kinetics and dynamics in water and aqueous solutions. Computer simulation study [53] first showed the non-exponential behaviour of the dynamics of forming and breaking H-bonds in liquid water due to the coupling between H-bond population and diffusion. Later, the jump reorientation mechanism [56] showed that water reorientation mechanism involves large amplitude angular jumps, rather than the commonly accepted sequence of small diffusive steps. Determination of H-bond lifetime is also a very important tool to examine the strength of H-bonds. This is largely affected by the structure of the local environment. It is important to mention here that very few attempts [61,62] have been made to calculate H-bond fluctuation dynamics for binary aqueous mixtures. Tetramethylurea (TMU) (shown in Scheme 1) is an amphiphilic molecule possessing both hydrophobic methyl groups and hydrophilic carbonyl group. TMU is easily soluble in water via H-bond interactions. TMU is an important biologically relevant molecule because it is a stronger denaturing agent than urea itself. [63,64] In this article, we show, through molecular dynamics simulations of TMU in water at various composition mixtures, that the confining effects of TMU are transduced on H-bond fluctuation dynamics of water. We have found that the longer lived water -water and water -TMU Hbonds exist in higher TMU concentrations. We have also observed that TMU induces significant influence on local environments that regulate the H-bond dynamics.
The organisation of the rest of the paper is as follows. Section 2 contains the molecular model and force field information. Simulation details are provided in Section 3. Simulation results and discussion are provided in Section 4. The paper ends with conclusions in Section 5.

Model and force field
Here we have considered the SPC/E model of water [65] and AMBER-type force field for TMU. The form of the force field is as follows: Here, U is the total potential energy of the system. The parameters b; b 0 ; u; u 0 ; f are actual bond length, equilibrium bond length, actual bond angle, equilibrium bond angle and dihedral angle, respectively. K b ; K u ; K f are the force constants for bond, angle and dihedral angle, respectively. fsc denotes the factor by which torsional barrier is divided. g is the phase shift angle and n the periodicity. The non-bonded interactions are included by the Lennard-Jones (LJ) potential. [66] The LJ parameters 1; s; r are the potential well depth, van der Waals radius and distance between atoms, respectively. The parameter q which represents the partial charge of the atom is taken to account for the Coulomb interaction. LJ interactions between unlike atoms are calculated by Lorentz -Berthelot combining rule, [66] where

Simulation details
Molecular dynamics simulations were performed using DL_POLY_Classic [67] package. We have used DL_FIELD [68] to make the input files in DL_POLY format. In this study, the methyl groups of TMU are treated by the united atom approach to avoid the complication of the simulation. All the force field parameters for TMU are taken from the literature. [69,70,71] A total of 512 molecules in a cubic box with periodic boundary condition have been considered. To obtain a reasonable initial configuration, we have run all the simulations in NPT ensemble at 1 atm pressure and 300 K using Hoover thermostat and barostat with relaxation times of 2 and 0.5 ps, respectively. We have run a total of 2 ns simulation in NPT ensemble to attain the experimental density of each system. Then, we have used the final configuration of NPT ensemble as the initial configuration for running simulations in NVT ensemble at 300 K using the Nosé -Hoover [72,73] thermostat with relaxation time of 0.5 ps. We have further equilibrated the system for a nanosecond to allow complete equilibration and then collected a 4 ns production run. Equations of motion were subsequently solved by leapfrog-Verlet algorithm [66] with time step of 1 fs. Verlet neighbour list shell width was set to 1.2 Å . Trajectories were saved after each 10 fs. Cut-off radius was set nearly to half of the box length. The Shake algorithm [74] was used to constrain all bonds involving hydrogen atoms. Electrostatic interactions were treated via the conventional Ewald summation technique. [66] In addition, we have used the TRANAL utility of MDynaMix program [75]  temperature are compared with those from experiments [76,77] in Figure 1. Good agreement between composition-dependent simulated and experimental densities indicates the validity of the force field employed here. We have also calculated 1 0 for water and TMU from simulated collective moment fluctuations and compared with the relevant experimental data. The following relation [78] has been employed to calculate 1 0 : with the collective moment, M ¼ P i m i , being the sum over individual dipole moments, m i , V the volume and k B T Boltzmann constant times the absolute temperature.
Our calculated 1 0 for water at 300 K temperature is , 73 and agrees well with experimental [79,80] and simulation [81] results. In contrast, the calculated 1 0 for TMU is found to be , 13 which is , 60% of the experimental value. [82] This discrepancy between simulation and experiment may arise partly from the united atom description for methyl groups in TMU and partly from system size. Note system size is an important factor in simulations for determining 1 0 as the latter is connected to the collective (k ! 0; k being the wavenumber) solvent modes. However, the H-bond dynamics that we intend to investigate is governed largely by the local environment and, therefore, consideration of 512 particles is expected to provide a semi-quantitative description of compositiondependent dynamics of this binary mixture.

Radial distribution functions
The relative structural arrangements of the atoms involved in H-bond interactions in the mixture are shown in Figure 2, where the simulated radial distribution functions (RDFs, gðrÞ) at several TMU mole fractions (X TMU ) are shown as a function of distance, r. While the upper panel depicts the simulated gðrÞ for oxygen and hydrogen atoms of water (O(W) -H(W)), the RDFs for water hydrogen and TMU oxygen atoms (O(TMU)-H(W)) are shown in the middle panel. At X TMU ¼ 0 (neat water), gðrÞ for O(W) -H (W) exhibits its first peak at , 1.7 Å , a distance nearly the sum of the radii of the atoms involved. Expectedly, the second peak appears at a distance (, 3.3 Å ) approximately equal to the sum of the diameters of these atoms. Note that the O(W) -H(W) gðrÞ peak intensity increases with the  increase of X TMU in the aqueous mixture, but the peak position remains nearly unaltered. This behaviour suggests that addition of TMU leads to enhanced interaction between oxygen and hydrogen atoms of water molecules. Interestingly, gðrÞ shown in the middle panel also indicate increased interaction between water hydrogen and TMU oxygen upon addition of TMU in water. Growth of both O(W) -H(W) and O(TMU) -H(W) RDFs with increase of TMU concentration suggests that water and TMU molecules arrange themselves in the binary mixture in such a way that the H-bonding interaction between water molecules and that between water and TMU molecules are facilitated upon addition of this amphiphile. This behaviour was observed earlier also for DMSO -water binary mixtures [44] and attributed to structural percolation.
It is interesting to note in these panels that the second peak intensity of O(TMU)-H(W) gðrÞ is much less than that of O(W) -H(W), suggesting more number of water molecules in the second solvation shell of water oxygen than that for TMU oxygen. The presence of methyl groups may be the reason for this relatively dwindled population of water molecules. This fact notwithstanding these RDFs does indicate successive enhancement of water -water and water -TMU interactions upon increase of X TMU in the aqueous mixture. This means weakening of interaction among TMU molecules in the presence of water. This is indeed the scenario as evidenced in the lower panel of Figure 2 which depicts H 3 C-CH 3 gðrÞ for TMU molecules in this mixture. Note in this panel the gðrÞ peak intensity decreases steadily and exhibits a trend opposite to that observed in the other two panels above. Importantly, this also indicates no aggregation among TMU molecules via methyl -methyl interaction even at large X TMU . However, the enhanced interactions between water -water and water -TMU suggest decrease of translational diffusion coefficients of these molecules with X TMU . This is what we have investigated next.

Mean squared displacements
The translational self-diffusion coefficient (D) of a particle moving in a fluid can be calculated from the mean squared displacements (MSDs) of centre-of-mass position vectors r i ðtÞ as follows: with kjDrðtÞj 2 l denoting the MSD. Subsequently, D has been calculated from the simulated kjDrðtÞj 2 l at long time as follows [83,84]: Simulated MSDs at various X TMU are presented in Figure 3 with upper panel depicting kjDrðtÞj 2 l for oxygen atom of water (upper panel) and lower panel for oxygen atom of TMU. Here we have presented MSD as a function of time in double logarithmic scale in order to show the t 2 and t 1 dependencies of kjDrðtÞj 2 l at short and long times, respectively. These dependencies are indicated in these plots by 'Slope I' and 'Slope II'. Note that there is no evidence of 'rattling in a cage' motion [85 -87] although simulated gðrÞ suggested stronger interactions between water and TMU molecules through H-bonding. Therefore, caging of TMU, even at low concentrations, by water molecules is probably non-existent in these aqueous binary mixtures.
Mixture composition dependence of viscosity (h) and translational diffusion (D) coefficients for water and TMU are investigated next. In simulations, h has been obtained by using the Green -Kubo relation, [83,84] h where a; b ¼ x; y; z. P ab denotes the off-diagonal terms ða -bÞ of the pressure tensor. Figure 4  comparison between the simulated and experimental [76] h for this mixture at , 300 K. Clearly, the agreement is semi-quantitative, indicating insufficiency of the model interaction potentials to represent quantitatively the interactions present in such aqueous binary mixtures, [81,88] and formation of aqueous complex through H-bond interaction. [89] The lower panel of Figure 4 shows the composition dependence of the simulated D for water and TMU for this mixture at 300 K. Simulated D for SPC/E water at , 300 K by other authors [81,88] is also shown to ensure fidelity of the present simulations. Note that D for both the specie first decreases with X TMU and then increases after a shallow minimum. Interestingly, one expects, from X TMU dependencies of medium density (see Figure 1) and viscosity, [76,89] a much sharper inverted parabolic dependence for D on composition for this mixture. Slip hydrodynamic predictions for D using both simulated (broken lines) and experimental (solid lines) viscosities are also shown in the same figure to highlight the difference between the simulations and hydrodynamic predictions. The difference between the simulated and predicted values may appear from the non-hydrodynamic character of the particle movement in these binary mixtures. Measurements of D for these mixtures are, therefore, necessary to ascertain whether the X TMU dependence of simulated D originates from partial decupling of diffusion from viscosity [90 -94] in these mixtures or the observed dependence is simply a reflection of the non-quantitativeness of the model interaction potentials used for accurate representation of the complex interactions in these aqueous binary mixtures.

Dynamic heterogeneity in (water 1 TMU) mixtures
Dynamic heterogeneity (DH) can impart non-hydrodynamic character to the particle movements, and can be understood in terms of non-Gaussian (NG) and new non-Gaussian (NNG) parameters, and single particle displacement distributions. [95 -97] The translational NG parameter (a 2 ) can be obtained from the simulated MSD by using the following relation: where kdr 2 ðtÞl ¼ kN 21 P N i¼1 jDrðtÞj 2 l and kdr 4 ðtÞl ¼ kN 21 P N i¼1 jDrðtÞj 4 l with dr denoting the single particle displacement. DrðtÞ ¼ r i ðtÞ 2 r i ð0Þ. Note that for homogeneous systems (for example, hot liquids), a 2 ðtÞ ¼ 0 for both at t ¼ 0 and t ¼ 1 but a 2 ðtÞ ¼ 0:2 at intermediate times. [90] Therefore, systems possessing DH would be characterised by a non-monotonic time dependence of a 2 ðtÞ with peak higher than 0.2. This peak signifies a DH timescale which is designated as t NG . Interestingly, another much slower DH timescale may also exist in the system which can be accessed via calculating the NNG parameter (gðtÞ) as follows [95]: with k1=dr 2 ðtÞl ¼ k1=N P N i¼1 1=jDrðtÞj 2 l. Clearly, gðtÞ is dominated by those particles which execute smaller displacements and contribute to a DH timescale (t NNG ) slower than that reflected by t NG . For systems possessing no spatially varying relaxation rates (that is, DH) t NG ¼ t NNG , and a 2 ðtÞ and gðtÞ curves would overlap.  In the lower panel, filled black circles represent simulated diffusion coefficients (D) for water and filled red triangles denote the simulated D for TMU. While the solid lines represent the slip hydrodynamic predictions by using the experimental h, the dashed lines denote the same but using the simulated h. For hydrodynamic calculations, radii of water and TMU molecules used are 1.425 and 3.5 Å , respectively. Star denotes the simulated diffusion coefficient for SPC/E water at 300 K reported in Ref. [79].
The DH signature can also be revealed by studying the single-particle displacement distribution as follows [95][96][97]: where G s ðdr; tÞ denotes the self-part of the van Hove correlation function. [83,84] Note for a Gaussian G s ðdr; tÞ, P½log 10 ðdrÞ; t becomes independent of time with a peak height of 2.13. Any deviation from the Gaussian shape and peak height would therefore indicate presence of DH in the system and may be interpreted as the origin for the partial decoupling of diffusion from solution viscosity. Figure 5 shows a 2 ðtÞ and gðtÞ simulated for water at four different mixture compositions at 300 K. Note that these compositions (that is X TMU ¼ 0.01, 0.2, 0.3 and 0.5) are chosen to represent dilute and concentrated aqueous solutions of TMU, including the composition that corresponds to the maximum in experimental and simulated h. Simulated a 2 ðtÞ and gðtÞ for neat water at 300 K, shown in Figure S1 of Supplementary information, depict a non-monotonic dependence with a peak value of , 0.2 at intermediate times, and the peak times (t NG and t NNG ) nearly overlap on each other. However, as X TMU increases, the peak value increases and t NG deviates from t NNG . The appearance of these two 'slow' timescales clearly demonstrates the presence of DH in the system. Therefore, partial decoupling between D and h for water noted in Figure 4 at higher TMU concentration may be attributed to the DH of these aqueous binary mixtures. Figure 6 depicts the simulated P½log 10 ðdrÞ; t for water at these representative mixture compositions at t NG and t NNG . Corresponding distributions for a Gaussian G s ðdr; tÞ are also shown for comparison. P½log 10 ðdrÞ; t at these DH timescales for ambient water is provided in Figure S2 (Supplementary information). Simulated distributions shown in these panels clearly demonstrate enhancement of DH for water upon increasing TMU concentration in these mixtures. This explains the partial decoupling of solution viscosity from water diffusion observed in Figure 4. Interestingly, the simulated P½log 10 ðdrÞ; t at t NG and t NNG for ambient water shows a slight deviation from the Gaussian character, which probably explains our simulated diffusion coefficient for neat water being , 10% larger than that reported in experiments (, 2.3 £ 10 25 cm 2 s 21 ).

H-bond dynamics
Simulated gðrÞ presented in Figure 2 has already revealed pronounced interaction between water and TMU in aqueous mixture, and this is expected to deeply impact the H-bond dynamics of water in solution. We have employed the following geometric definition of H-bonding between a pair of molecules [57 -60]: 1. Distance between two oxygen atoms must be less than 3.5 Å .

Angle between OZO bond vector and OZH bond
vector must be less than 308. 3. Distance between O and H must be less than 2.45 Å .
How long a pair of molecules remains H-bonded can be calculated by two different time-correlation functions. [23 -26] The relaxation of H-bond lifetime is calculated by the time-correlation function, S HB ðtÞ, where HðtÞ ¼ 1 if the tagged pair of molecules, for which hð0Þ is calculated, remains continuously H-bonded till time t, else zero. S HB ðtÞ describes the probability that a tagged pair of molecules remains continuously H-bonded up to time t, and approaches to zero when H-bonding between them breaks down for the first time. The time profile of S HB ðtÞ therefore represents relaxation of H-bond lifetime, and hence, the average H-bond lifetime, kt HB S l, can be obtained via time integration of this normalised time correlation function.
The structural relaxation time of H-bond is calculated by the time-correlation function, C HB ðtÞ; defined as [23 -26] where C HB ðtÞ exhibits a much slower decay than S HB ðtÞ; as the former is independent of possible breaking of H-bonds between the tagged pair of molecules in the intermediate time and considers the reformation of bond after breakage. C HB ðtÞ is, therefore, associated with H-bond structural relaxation, and thus the time integration of this normalised correlation function provides the average structural relaxation time, kt HB C l. Figure 7 shows the decay of S HB ðtÞ and C HB ðtÞ as a function of time simulated at X TMU ¼ 0 and T ¼ 300 K. Bi-exponential fit to S HB ðtÞ and tri-exponential fit to C HB ðtÞ are also shown in these figures. Note simulated S HB ðtÞ and C HB ðtÞ at other mixture compositions can also be similarly described via bi-and tri-exponentials, respectively. Fit parameters required for these correlation functions at all X TMU values are summarised in Tables S1 -S4 of the Supplementary information. Note that multi-exponential X TMU =0.01 nature of S HB ðtÞ and C HB ðtÞ has been found earlier also for a variety of systems, ranging from aqueous micellar solutions [60,98,99] and cyclodextrins [100] to moisturised proteins and DNA. [101,102] Note kt HB S l and kt HB C l for pure water are, respectively, , 0.5 and , 4 ps, which are in good agreement with earlier findings. [59 -61,98 -102] The multi-exponential character of the relaxations may have origin in the two different density regions that exist in liquid water -one is low density four-coordinated tetrahedral H-bonding region and the other is two or threecoordinated high density region. [103,104] However, further study is required to determine whether density fluctuations in these differently populated regions are indeed responsible for the multi-exponential relaxations of S HB ðtÞ and C HB ðtÞ in neat water and in water-rich solutions. Figure 8 presents the X TMU -dependent relaxations of S HB ðtÞ for water -water and water -TMU H-bond lifetimes. Note S HB ðtÞ for both cases slows down with each successive addition of TMU in the mixture, resulting in longer H-bond lifetimes in the presence of TMU. Interestingly, the enhancement of water -water H-bond lifetime upon gradual addition of TMU in aqueous solution resembles the slowing down of water confined in carbon nanotube [105] and water at the surface of an antifreeze protein. [106] Therefore, it may be stated that the presence of TMU slows down the H-bond dynamics in some ways similar to those operative inside the confinement and at the micelle and protein surfaces. One of the frames from the simulation trajectory at X TMU ¼ 0.2 is shown in Figure 9 which depicts the alignment of water and TMU molecules in the binary mixture. Figure 9 indeed suggests that water molecules near the TMU surface can mimic the effects observed either under confinement or near a micellar surface. It is to be mentioned here that longer H-bond lifetime does not necessarily indicate stronger H-bonds. [107] The average lifetime of H-bonds largely depends on local environments. From the simulated gðrÞ we have noticed growth of O -H interactions between water -water and water -TMU molecules. This structural compactness probably forces the water and TMU molecules to remain H-bonded for longer times. Quantum mechanical calculations of H-bond energies may provide further information on this aspect.
H-bond structural relaxation functions, C HB ðtÞ, are shown in Figure 10 for water -water (upper panel) and water -TMU (lower panel) specie. As found for S HB ðtÞ in Figure 8, C HB ðtÞ also slows down upon addition of TMU in the aqueous solution, producing longer average structural relaxation times, kt HB C l. This can also be linked to the changes in gðrÞ upon addition of TMU in the mixture. The  X TMU -induced slowing down of S HB ðtÞ and C HB ðtÞ relaxations is summarised in Figure 11 where average relaxation times, kt HB S l and kt HB C l, are shown as a function of TMU mole fraction in the mixture. Note that both these relaxation times increase linearly with X TMU , and do not show any non-monotonic composition dependence. This clearly indicates that none of these relaxation processes are governed by medium viscosity. Moreover, kt HB S l and kt HB C l for water -TMU are longer at all X TMU than those for water -water. Relatively lower diffusion of TMU may not be the reason for this as these relaxation times do not couple to solution viscosity. Data in Figure 11 also indicate that kt HB S l and kt HB C l slow down by a factor of , 3 -10 upon changing the medium from neat water to aqueous TMU solution with X TMU ¼ 0.75. Interestingly, this extent of slowing down is comparable to lengthening of these quantities for water encapsulated in cyclodextrin cavities. [100]

Conclusion
In summary, this simulation study suggests considerable slowing down of H-bond relaxation dynamics in the presence of amphiphilic molecules in bulk binary mixtures, and the extent of slowing down mimics those found for encapsulated water molecules in cavities and those near micellar and biological surfaces. Simulated RDFs indicate increase of water -water and water -TMU interactions upon addition of TMU in the aqueous solution, but decrease of TMU -TMU interaction. Simulated MSDs do not show any signature of cagerattling and this has been interpreted as no evidence for the formation of TMU -water clathrate-type structure even at low TMU concentration. Simulated diffusion coefficients for both the specie seem to indicate partial decoupling of diffusion from medium viscosity. Simulated DH and single particle displacement distributions reflect deviations of particle movements from hydrodynamic behaviour. In addition, simulated diffusion coefficients do not show as pronounced composition dependence as that shown by viscosity coefficients. Use of polarisable potential and explicit treatment of methyl group in the force field for TMU may improve the agreement between the composition-dependent simulated and experimental viscosity coefficients. Relaxation profiles of H-bond lifetime and structure are multi-exponential in character, and we propose that density inhomogeneity inherent to liquid water may be responsible for this non-exponential behaviour. This aspect should be further investigated. In addition, DH of water molecules, particularly of those   Figure 11. (Colour online) Plot of average relaxation times for X TMU -dependent H-bond dynamics and structural dynamics. As discussed in the text, kt HB S l has been obtained analytically from fits of S HB ðtÞ and kt HB C l similarly from C HB ðtÞ. Solid circles in both panels represent simulated data reported in Ref. [59]. present at high TMU concentration, needs to be explored as the simulated diffusion coefficients appear to suggest decoupling from medium viscosity.