How do glycerol and dimethyl sulphoxide affect local tetrahedral structure of wateraround a nonpolar solute at low temperature? Importance of preferential interaction
Snehasis Daschakraborty

Citation: The Journal of Chemical Physics 148, 134501 (2018); doi: 10.1063/1.5019239
View online: https://doi.org/10.1063/1.5019239
View Table of Contents: http://aip.scitation.org/toc/jcp/148/13
Published by the American Institute of Physics

THE JOURNAL OF CHEMICAL PHYSICS 148, 134501 (2018)

How do glycerol and dimethyl sulphoxide affect local tetrahedral
structure of water around a nonpolar solute at low temperature?
Importance of preferential interaction
Snehasis Daschakrabortya)Department of Chemistry, Indian Institute of Technology Patna, Patna-801103, Bihar, India(Received 12 December 2017; accepted 12 March 2018; published online 2 April 2018)
Glycerol and dimethyl sulphoxide (DMSO) have vital roles in cryoprotection of living cells, tissues,
etc. The above action has been directly linked with disruption of hydrogen (H-) bond structure anddynamics of water by these cosolvents at bulk region and around various complex units, such as
peptide, amino acid, protein, and lipid membrane. However, the disruption of the local structure of
the water solvent around a purely hydrophobic solute is still not studied extensively. The latter is also
important in the context of stabilization of protein from cold denaturation. Through all-atom molecular
dynamics simulation, we have investigated the comparative effect of glycerol and DMSO on theorientational order of water around a nonpolar solute at 5 ◦ C. A steady reduction of the tetrahedral
order of water is observed at bulk (>10 Å distance from the solute) and solute interface (<5.5 Å
distance from the solute) with increasing the cosolvent concentration. Contrasting roles of glycerol
and DMSO have been evidenced. While DMSO affects the H-bond structure of the interfacial water
more than that of the bulk water, glycerol affects the water structure almost uniformly at all regions
around the solute. Furthermore, while glycerol helps to retain water molecules at the interface, DMSO
significantly reduces the water content in that region. We have put forward a plausible mechanism
for these contrasting roles of these cosolvents. The solute-cosolvent hydrophobic-interaction-inducedorientational alignment of an interfacial cosolvent molecule determines whether the involvement of
the cosolvent molecules in H-bonding with solvent water in the interface is akin to the bulk region ornot. Published by AIP Publishing. https://doi.org/10.1063/1.5019239

I. INTRODUCTION

Cosolvents, like glycerol and dimethyl sulphoxide(DMSO), in liquid water protect living cells and tissues by
preventing freezing of water at extremely low temperature.
Thus the cosolvents are popularly called as cryoprotectants,
and the process of prevention from freezing is known as cryopreservation. Several experimental and computer simulation
studies indicated that the anti-freezing action of these cryoprotectants stems from their abilities in disrupting the threedimensional tetrahedral H-bonding network and underlying
molecular dynamics (MD) of liquid water.1–21 It was observedthat the hydroxyl groups of glycerol molecules form strong
H-bonds with water and the hydrophobic residues aggregate
together. Thereby the regular three-dimensional H-bond network of water is severely disturbed. However, the scenario
can be astonishingly different at molecular length scale. For
example, it was reported that water forms lower density, moreordered tetrahedral water structure in high glycerol concentration than pure water at room temperature.2–4 This might
be due to the nano-segregation of glycerol-rich and waterrich regions in the mixture, which protects water to form alow-density by an extensive and encapsulating glycerol interface. Multiple experimental and simulation reports indicateda)Author to whom correspondence should be addressed: snehasis@iitp.ac.in

0021-9606/2018/148(13)/134501/11/$30.00

similar disruption of three-dimensional H-bonding of water inthe presence of DMSO.5–10,15,16,19–21
Notwithstanding the fact that the above studies provide a
broad view on structural and dynamical perturbation of waterby these cosolvents at the bulk regime, they offer limited
insight on how these cryoprotectants affect water H-bonding
structure and dynamics near a purely hydrophobic solute. The
latter is particularly important in understanding the mechanismof protection of proteins by these cosolvents from cold denaturation. The cold denaturation is often linked with the weakhydrophobic interaction among the hydrophobic residues
of protein at low temperature.22–49 However, the increased
H-bond interaction between polar residues of amino acid andsolvent water at low temperature may also play an important role in cold denaturation of protein.50–52 Apart from theabove two pictures, another exciting view also exists. A Monte
Carlo simulation-based approach indicated that the increaseof water· · · water H-bonding at the protein interface with lowering of temperature is the sufficient condition for protein’s
cold denaturation.53 We will return to this discussion again in
Sec. III.
Some reports are found in the literature on how glycerol
and DMSO affect the water solvent structure and dynamics
near peptide, amino acid residues, protein, lipid membrane,
etc. It was found that both the cosolvents substantially perturb
water solvent structure and dynamics near all the above entities. However, the mechanisms of the above actions of these

148, 134501-1

Published by AIP Publishing.

134501-2

Snehasis Daschakraborty

two cosolvents are notably different. The MD simulation study
by Johnson et al. has analyzed the role of the two cosolvents in
peptide hydration structures.54 While the presence of glycerol
in the solution increases hydration number near the peptide
interface, DMSO effectively reduces it. Similar contrasting
behavior was seen on the lipid membrane interface.55 Glycerol influences the diffusivity of water almost uniformly near
the bilayer surface and that in the bulk solution, but DMSO
increases the diffusivity of surface water more than that of
the bulk water molecules. It is also known that glycerol and
DMSO act very differently for stabilization of protein. Glycerol is an excellent protein stabilizer,54,56,57 whereas DMSO’s
action is diverse and strongly dependent on the concentrationof DMSO in the solution.58–64 For example, with increasing
DMSO concentration, initially (at x D < 0.05) the conformational flexibility of lysozyme increases, followed by the slight
unfolding of native state in the DMSO concentration range
0.10 < x D < 0.15, and finally unfolding of the protein afterx D > 0.15.58 This picture was later supported by an MD
simulation study.62
Although the above studies provide significant insights onsome specific questions, a more fundamental question remains
mostly unaddressed. The question is, how do these cosolventsaffect the strong H-bonding clathrate like cage structure ofwater 65,66 —popularly known as “iceberg,” first proposed byFrank and Evans and later evidenced by several experimental and computer simulation studies67–73 —around a purelyhydrophobic solute?
We have studied here the effect of glycerol and DMSO
concentration on the orientational ordering of solvent watermolecules around a very simple nonpolar solute, xenon (Xe)using classical MD simulation. It is well known that Xe is
prototypical hydrophobic solute and its hydration structureand dynamics in liquid water resemble other popularly knownhydrophobic solutes, like methane and ethane.67,73–77 However, one can still ask whether the influence of the cosolventson the hydration structure around the small nonpolar solute,Xe, is analogous to the influence on the hydration structure around much more complicated and bigger entities, suchas protein and lipid membrane. It was observed earlier thatthe hydration structure around a hydrophobic solute changes
significantly with an increase in the size of the solute.78–84
The pioneering Lum-Chandler-Weeks theory predicts that forsmall solutes the free energy of hydration of a hard-spheresolute increases almost linearly with the solute volume but
beyond a critical length scale (also termed as small-to-large
crossover length scale) the free energy increases linearly withthe surface area of the solute.78 The theory, therefore, predicts drying of the extended nonpolar surface only beyond thesmall-to-large crossover length scale (∼1 nm) of the nonpolar
solute.
A molecular simulation was also carried out to study theeffect of various thermodynamic properties and concentrationof additive, like ethanol, on the above small-to-large crossoverlength scale.81 (Note that this length scale is ∼1 nm understandard temperature and pressure condition.) It was foundthat the crossover length scale can decrease up to the size of asmall molecule after adding ethanol like cosolvent. This is very
important in the context of the present study because we have

J. Chem. Phys. 148, 134501 (2018)

also used the cosolvents, glycerol or DMSO, which are notvery different from ethanol regarding interaction with water.
Therefore we expect a decrease of the small-to-large crossover
distance to the size of the small solute, like Xe, for the mixtures. Our choice of the solute in the present study is further
emboldened by a recent molecular simulation study,85 which
intended to elucidate whether proteins behave effectively likesmall or large non-polar solute (in terms of the hydrophobic
effect). The simulation has reported the following observation.
Despite the proteins are large units (larger than 1 nm), they
can behave effectively like “small” particles (i.e., smaller than
1 nm) in terms of hydrophobic effect. The complex structure ofproteins, having both polar and non-polar residues, is the possible reason for the above behavior. Therefore, the proximityof polar and nonpolar residues in protein makes protein hydration structure very similar to small nonpolar solute, which
influences the solvent water H-bonding strength only slightly.
Therefore, we expect that the observation for the system canbe transferred to systems of practical interest for understanding the role of cryoprotectants in influencing the hydrophobic
structure.
The simulation was done at 268 K temperature given theinterest in the cosolvent activity at that temperature aroundwhat cold denaturation often occurs for various proteins. In
this study, we are asking the following questions: How dothese two cosolvents perturb the local tetrahedral structureof water around the solute? Do the cosolvents uniformly
affect the structure of water inside and outside the first hydration shell? Are there any contrasting effects of glycerol and
DMSO on hydrophobic hydration structure? Is there any roleof preferential solvation?
The outline for the remainder of this paper is as follows. The molecular models and other simulation details aredetailed in Sec. II. Section III presents the simulation results
and discussions. Concluding remarks are offered in Sec. IV.
II. SIMULATION DETAILSIn order to study the composition dependence of theresults, we have considered here 6 compositions (with themole fraction of glycerol or DMSO, x = 0.02, 0.05, 0.07,
0.10, 0.12, and 0.15) for each of the glycerol+water (GW) andDMSO+water (DW) mixtures plus neat liquid water, resulting
in total 13 different systems. We have not gone beyond themole fraction of the cosolvent x = 0.15 because of the following reason. Although protein is stable at a far wide range ofcompositions of GW mixtures (including x G = 0.99!),86 theformer is stable only at a narrow range of composition in the
DW mixture. It was observed that lysozyme starts denaturing
beyond x D = 0.15 mole fraction of DMSO in the DW mixture.58–64,87 Therefore, we have limited our studies up to thecomposition, x = 0.15, where the protein stays stable in its
native state.
All-atom modeling of glycerol molecules are done using
AMBER88 type force field parameters, initially developedby Chelli et al.89 and later modified by Blieck et al.90 The
TIP4P/2005 model91 has been used for water molecules. This
combination of the force field models is proven to be themost successful in reproducing experimental thermodynamic

134501-3

Snehasis Daschakraborty

properties and other microscopic structural and transport properties of the GW binary mixture at a wide range of temperatures including sub-zero.92 We have used here the fully
atomistic Strader-Feller model of DMSO,93 which is highly
successful in predicting experimentally observed thermodynamic and microscopic quantities. The nonpolar solute Xe ismodeled as 12-6 Lennard-Jones (LJ) particle having LJ parameters used elsewhere.74,75 The solute–solvent cross parameters are deduced from the Lorentz–Berthelot mixing rules,
which is the best combination rule for the aqueous solutionof Xe.76
First, we optimized a glycerol and a DMSO molecule separately in the gas phase at the B3LYP/6-311+G∗∗ level oftheory94,95 and basis set using the Gaussian 09 package.96
Then we prepared 13 simulation boxes—each having box
length of 30 Å—by mixing an appropriate number of geometrically optimized glycerol and DMSO molecules separately
with water following the composition of the solvent mixture
and one Xe solute. Each simulation box contains 500 solvent
molecules and 1 solute molecule. The initial configurations of
the molecules were set up randomly using the PACKMOL software.97 The compositions of the mixtures have been detailedin Table S1 of the supplementary material.
MD simulations for the above 13 systems were then initiated using the Gromacs software package.98 All the systems
were first equilibrated using the NPT ensemble for 10 ns.
The temperature was maintained at 268 K using the NośeHoover thermostat (coupling time 0.5 ps)99,100 and pressure at
1 bar using the Berendsen barostat (coupling time 0.5 ps).101
In order to confirm that the Berendsen barostat is working correctly—specifically, whether the resulting density ofthe system remains same on using different pressure coupling methods102 —we re-equilibrate the systems with theParrinello-Rahman barostat 103 with the same coupling constant of 0.5 ps. The equations of motions were always solved
using the Verlet leapfrog algorithm every after 1 fs. Table S1of the supplementary material lists the final simulated densities separately for the above two barostats. Clearly, the final
densities are almost the same, indicating consistent density
equilibration with the two barostats. Also, the final densities
are in excellent agreement with earlier simulation and experimentally measured densities.92,104 This further supports our
choice of force field models for the solvent molecules and othersimulation parameters.
For obtaining production trajectories, NVT simulations
were initiated from the final configurations of the previous
NPT simulations. These simulations were run for 40 ns bysolving equations of motion every after 1 fs using the Verlet
leapfrog algorithm. The temperature constraining was done byusing the Nośe-Hoover thermostat, one of the efficient temperature control algorithms,102 at 268 K with a coupling constantof 0.5 ps. The configurations were saved after every 50 fs time
for analysis. The above simulation protocols are akin to otherrelevant simulation studies.
III. RESULTS AND DISCUSSIONSIn this section, we first discuss the consequence ofincreasing glycerol and DMSO concentration in aqueous

J. Chem. Phys. 148, 134501 (2018)

solution on the local tetrahedral order of liquid water aroundthe hydrophobic solute. We continue our discussion in
Subsections III B–III C about the roles of these two cosolvents
in H-bonding and spatial structure of water solvent around thesolute. Finally, we illustrate the observations by looking at theorientational arrangement of the cosolvent molecules in thehydration shell of the Xe solute.
A. Tetrahedral order parameter

Orientational tetrahedral order parameter (q) is the mostextensively used parameter for understanding the tetrahedral
structure of a liquid. (See, e.g., Refs. 67 and 105–110.) For neatliquid water, q is calculated by considering 4 nearest water
oxygen (OW ) atoms from a central water OW and using thefollowing equation:!2
1
3 X3 X4.
(1)
cos θ jk +
q =1−
j=1
k=j+1
8
3
Here, θ jk is the angle formed by the lines joining OW of agiven water molecule and those of its four nearest neighbors j
and k. For entirely random configuration, like an ideal gas, theaverage value hqi is 0. On the other hand, hqi = 1 for perfectly
tetrahedral configuration, like ice.
First, we study q for water clusters in the two mixtures.
Figures 1(a) and 1(c) present the normalized probability density distributions of q for pure water, 3 GW mixtures (x G = 0.05,
0.10, and 0.15), and 3 DW mixtures (x D = 0.05, 0.10, and 0.15).
In this calculation, we have considered only the nearest fourwater molecules from a central water molecule. The rationale
is to check how the cosolvents prevent tetrahedral orderingof water among themselves around the hydrophobic solute at
low temperature. q is calculated for 1st hydration shell (interfacial) and bulk water molecules separately. The interface andthe bulk are identified from radial distribution function (RDF)
between the solute Xe and water OW , presented in Figs. S1and S2 of the supplementary material. The interfacial layer0extends up to 5.5 Å distance from the solute (position of the0
first minima of Xe−−OW RDF). The space outside 10 Å distance from the solute (Xe−−OW RDF almost levels off after
this distance) is considered as bulk. Note that these distancecriteria are consistent with a recent MD simulation study.111
Figures 1(a) and 1(c) show that the distribution of q—both atinterface and bulk—is maximum at q ∼ 0.83 with a shoulder
at q ∼ 0.50 for pure water. The local tetrahedral structure ofwater at the present simulation temperature (T = 268 K) is
palpably more ordered than the structure at T = 298 K.67,110
The intensity of the peak for both the interfacial and bulk
water decreases steadily. (As suggested in Ref. 53, this lowtemperature ordering of solvent water near nonpolar solute
can be the key for cold denaturation of the protein.) However,no significant shift of the peak is observed. Closer inspectionof Figs. 1(a) and 1(c) further reveals that while the intensityof the peak—centered at q ∼ 0.83—gradually decreases, newpeak burgeons at q ∼ 0.45. The intensity of the new peakprevails over the intensity of the peak at q ∼ 0.83 beyond
x = 0.10.
In order to elucidate further, we view the structures of
5 water clusters including the central water for q ∼ 0.83 and134501-4

Snehasis Daschakraborty

J. Chem. Phys. 148, 134501 (2018)

FIG. 1. Normalized probability distribution of the tetrahedral order parameter, q, for bulk (dashed lines) and interfacial water molecules (solid lines) forneat water and binary mixtures of waterand cosolvents at different compositions: (a) GW mixtures and (c) DW mixtures. Difference between the distributions of q in bulk and the first hydration
shell around the solute for (b) GW and(d) DW binary mixtures.

q ∼ 0.45 for both GW and DW mixtures. Four such water clusters are shown in Fig. 2. Apparently, the structures in Figs. 2(a)and 2(b) are close to the regular tetrahedron. Therefore, thepeak at q ∼ 0.83 corresponds to a tetrahedral-like arrangementwith the four water molecules situated at the four vertices of
the tetrahedron. On the other hand, Figs. 2(c) and 2(d) revealthat the peak at q ∼ 0.45 corresponds to a broken tetrahedral geometry with one water molecule dislocated from one
vertex of the tetrahedron. This essentially suggests that thegradual burgeoning peak at q ∼ 0.45 with increasing glycerol/DMSO concentration—in Figs. 1(a) and 1(c)—is due tothe disruption of the local tetrahedral geometry of water atboth the bulk and hydrophobic interface. At x > 0.10, the water
clusters turn out to be dominated by the distorted tetrahedral
structure. We will show later that the fourth free vertex of
the tetrahedron is actually occupied by a cosolvent moleculeand thereby uphold “hybrid” tetrahedral clusters (composed ofboth water and cosolvent molecules) around the center watermolecule.

Now, we turn our focus on comparative analysis between
the distributions of q for interfacial and bulk water molecules.
The differences between the two distributions are shown in
Figs. 1(b) and 1(d) for GW and DW mixtures, respectively.
For pure water, the tetrahedral structure of interfacial water
solvent is more ordered than that of the bulk. This is consistent
with the picture that the strong H-bonding structure of interfacial water exists in the vicinity of the hydrophobic solute.67–73
Figures 1(a) and 1(c) demonstrate that even though the tetrahedral orders of both the bulk and interfacial water molecules
decrease with increasing mole fraction of glycerol in the GW
mixtures, the interfacial water molecules always stay moreordered than bulk water. However, the picture is markedly
different for DW mixtures. The interfacial water molecules
become eventually less ordered than the bulk waters at
x D > 0.05 for DW mixtures.
To elucidate further with the more quantitative analysis,
we plot average tetrahedral order parameters of the interfacial
water molecules (hqI i) and those of the bulk water molecules

FIG. 2. Representative water clusters for q ∼ 0.83 andq ∼ 0.45 for GW and DW binary mixtures.

134501-5

Snehasis Daschakraborty

FIG. 3. Average tetrahedral order parameters of the interfacial watermolecules (hqI i) (solid line with filled circles) and those of the bulk watermolecules (hqB i) (dashed line with open circles) in the two mixtures: (a) GW
and (b) DW as a function of x G and x D , respectively. (c) The ratios hqI i/hqB i
as functions of x G and x D for the two mixtures. The red line with open circles and the blue line with filled circles represent for GW and DW mixtures,
respectively, in panel (c). Horizontal dashed blue line in panel (c) indicateshqI i/hqB i = 1.

(hqB i) for GW and DW mixtures in Figs. 3(a) and 3(b), respectively. Evidently, both hqI i and hqB i decreases almost linearly with increasing the cosolvent concentration. The ratio
hqI i/hqB i is plotted in Fig. 3(c) as functions of x G and x D
for GW and DW mixtures, respectively. The ratio hqI i/hqB i
is above unity for all the compositions of the GW mixturesand stays almost insensitive on glycerol concentration up tox G = 0.15, the maximum glycerol concentration in the present
study. This immediately suggests that glycerol uniformly
affects the interfacial and bulk water. However, the scenario is
entirely different for DW mixtures. Figure 3(c) shows that hqI i
decreases more rapidly than hqB i with increasing x D in DW
mixtures. Therefore the ratio hqI i/hqB i becomes less than unity
at x D > 0.05 and decreases up to a value of 0.9 at x D = 0.15. Thissuggests that, unlike glycerol, DMSO affects the tetrahedral
arrangement of the interfacial water more radically than that of

J. Chem. Phys. 148, 134501 (2018)

the bulk water. The observation that glycerol uniformly influences the tetrahedral structure of water around the hydrophobic
solute and that DMSO disrupts the tetrahedral arrangement ofinterfacial water more drastically than that of the bulk waterstrongly recommends the contrasting roles of these two cosolvents in affecting the tetrahedral structure of water arounda hydrophobic solute. Similar contrasting behaviors of thesetwo cosolvents were reported earlier while studying the effectof these cosolvents on hydration structure and dynamics near
the peptide water interface.54 While glycerol increases thehydration number near the peptide interface, DMSO decreases
the hydration number near the same peptide interface. Thismight be due to the preferential hydration of the hydrophobic residues. We will return to this discussion in Sec. III B. A
more recent experimental study has also revealed that DMSO
and glycerol affect the diffusivity of interfacial and bulk water(from lipid bilayer surface) very differently.55 While glycerol
uniformly affects the diffusivity of water irrespective of the
distance from the lipid bilayer surface, DMSO significantly
increases the diffusivity of surface water compared to that of
the bulk water.
Now, we turn our focus on how glycerol and DMSO
cosolvents directly participate in the local tetrahedral structure of liquid water. While glycerol molecules can act as boththe H-bond donor and acceptor with water solvent molecules,
DMSO molecules only accept H-bonds from the solvent water.
In the beginning of this section, we have indicated that one
vertex—in majority of the tetrahedral clusters of water aroundthe hydrophobic solute—turns out to be liberated from water
when the mole fraction of the cosolvent reaches x = 0.15 andthat a cosolvent molecule occupies this unattended position byforming H-bond with the central water solvent molecule. Figures 4(a) and 4(b) present the normalized probability density
distributions of q for pure water, 3 GW mixtures (x G = 0.05,
0.10, and 0.15), and 3 DW mixtures (x D = 0.05, 0.10, and0.15). In this calculation, we have considered the nearest foursolvent molecules (both water and cosolvents) from a central
water molecule. Clearly, the distributions for both GW and DW
mixtures culminate at q ∼ 0.8 with a shoulder at q ∼ 0.5. Unlike
in Fig. 1, the peak of the distribution remains fixed at q ∼ 0.8
for all compositions and no new peak emerges at q ∼ 0.5 except
there is a slight rise of the shoulder. Clearly, the peak intensity
decreases with increasing cosolvent concentration. However,the loss of the peak intensity is much less than when only water
clusters are considered in the calculation. This indicates thatthe retention of the tetrahedral structure occurs at least up tox = 0.15. However, the pure water tetrahedral clusters turn into
the multicomponent tetrahedral cluster. These “hybrid” tetrahedral structures may potentially reduce the probability of ice
crystal formation and their growth. These observations support the fact that these cosolvents act as potential anti-freezing
substances.1–21
We have presented, in Figs. 5(a) and 5(b), the averagetetrahedral order parameters hqI i and hqB i for the interfacial
and bulk water solvent molecules, respectively, in DW andGW mixtures. These are calculated for the tetrahedral clusters
composed of both water and cosolvents. Figure 5(a) shows thatboth the hqI i and hqB i decrease almost parallel with increasing x G . Naturally, the ratio hqI i/hqB i remains almost constant134501-6

Snehasis Daschakraborty

J. Chem. Phys. 148, 134501 (2018)

FIG. 4. Distribution of the tetrahedral order parameter, q, for bulk (dashedlines) and interfacial water molecules (solid lines) for neat water and binary
mixtures of water and cosolvents at different compositions: (a) GW mixturesand (b) DW mixtures. Note that these calculations are done by considering
the nearest four solvent molecules (both water and cosolvents) from a central
water molecule.

at slightly higher than unity. The scenario is again different
for the DW mixture. As visible in both Figs. 4(b) and 5(b),hqI i at the interface decrease more than hqB i in bulk with thesame increase of x D . Therefore the ratio, presented in Fig. 5(c),decreases from slightly higher than 1.0 for neat water to ∼0.95
at x D = 0.15. The above observations clearly indicate that while
glycerol affects the tetrahedral order of the interfacial and thebulk water solvent almost uniformly, DMSO’s influence is not
uniform. A more meticulous comparison between Figs. 5(a)and 5(b) reveals an interesting picture. At x G = 0.15, the valueof hqI i at the interface (0.62) is higher than the value (0.59) atthe same mole fraction for the DW mixture. On the other hand,
the hqB i value (0.62) for GW mixture at x G = 0.15 is almost
the same with the hqB i value for DW mixture at x D = 0.15.
Consequently, the glycerol and DMSO cosolvents perturb thethree-dimensional tetrahedral H-bonding network of water
near the hydrophobic solute almost similarly at the bulk limit,
but the effect on the interfacial water solvent is not akin. DMSO
disrupts the orientational order of interfacial water more
intensely. This observation gives the signature of stronger possible preferential interaction between the hydrophobic soluteand DMSO cosolvent. We will return to this discussion in
Sec. III D.
B. H-bonding of water

In order to elucidate further on how the two cosolvents
relatively affect H-bonding between water molecules at thehydrophobic interface and the bulk, we calculate the numberof H-bonds per water molecule (N h ) as a function of distancefrom the solute.

FIG. 5. Average tetrahedral order parameters of the interfacial watermolecules (hqI i) (solid line with filled circles) and of the bulk water molecules(hqB i) (dashed line with open circles) in the two mixtures: (a) GW and (b) DW
as a function of x G and x D , respectively. (c) The ratios hqI i/hqB i as functionsof x G and x D for the two mixtures. The red line with open circles and theblue line with filled circles represent for GW and DW mixtures, respectively.
Horizontal dashed blue line in panel (c) indicates hqI i/hqB i = 1. Note thatthese calculations are done by considering the nearest four solvent molecules(both water and cosolvents) from a central water molecule.

Figures 6(a) and 6(b) exhibit N h as a function of distancefrom the solute for 4 different compositions (x = 0.0, 0.05, 0.10,
and 0.15) for each of the GW and the DW binary mixtures.
We have adopted the widely accepted H-bond criteria (rOD OA< 3.5 Å, rOA HD < 2.45 Å, HD OD OA < 30◦ , where “D” and“A” stand for the H-bond donor and acceptor, respectively).
Noticeably, for neat liquid water, N h goes to the maximum atrXeOW ∼ 4.5 Å distance from the solute and then after someundulation it almost levels off beyond rXeOW = 10 Å distance.
The peak at rXeOW ∼ 4.5 Å distance indicates stronger water
H-bonding network at the solute-water interface. This result
is consistent with the earlier observations of stronger water
H-bonding near the hydrophobic interface than that in thebulk.73 Figure 6(a) shows that for GW binary mixtures, N hdecreases steadily with increasing x G at all the distances fromthe solute. To see the relative perturbation by the cosolventsat the interface and the bulk regions, we have plotted theratio NhI /NhB in Fig. 6(c). Here, NhI and NhB are, respectively,134501-7

Snehasis Daschakraborty

J. Chem. Phys. 148, 134501 (2018)

decreases with an increase in x D . This indicates that DMSO
molecules disturb the H-bond structure of the interfacial watermolecules more strongly than that of the bulk water moleculesin the DW mixture. Also, as we have seen in Fig. 3, the local
tetrahedral structure of hydrophobic interfacial water is disrupted more than that in bulk. Therefore, DMSO decreases
both the number of water· · · water H-bonding and the local
tetrahedral H-bonding network. This can be one of the possible reasons of why protein gets precipitated in high DMSO
concentration in solution.
We have repeated the above analysis but considering bothwater and cosolvent as H-bond partners of the central watermolecule. Similar results are obtained. The results are presented in Fig. S3 of the supplementary material. It is worth
mentioning that the decrease of the overall H-bonding interaction with increasing the cosolvent concentration can explain
the stability of hydrophobic-interaction-induced self-assembly
structures, like proteins. This is due to the theory, describedin Ref. 53, which states that cold denaturation of protein canoccur naturally because of increasing strength of interfacial
water· · · water hydrogen bond interaction.
C. Radial distribution of water solvent aroundthe solute

In Secs. III A and III B, we found strong indication thatglycerol draws more water molecules at the solute/water interface from bulk and DMSO excludes water solvent from theinterface. To obtain more conclusive evidence, we look at themass density of water (ρw ) as a function of distance fromthe solute for GW and DW mixtures. Figures 7(a) and 7(b)present ρw as a function of distance from the solute at 4 different compositions of GW and DW mixtures, respectively.
FIG. 6. [(a) and (b)] Number of water· · · water H-bonds per water molecule(N h ) as a function of distance from the solute for GW mixtures [panel (a)]and DW mixtures [panel (b)]. The data for only 4 different compositions(x = 0.0, 0.05, 0.10, 0.15) are shown here. (c) The ratio between N h at theinterface (N h I ) and N h at the bulk (N h B ). The solid lines in panel (c) are
the linear fits of the data. The vertical lines in panels (a) and (b) indicate theupper limit of the interface (1st hydration shell) and the lower limit of thebulk region. The horizontal dashed lines in panels (a) and (b) indicate the bulk
values of the corresponding compositions of the two mixtures. The horizontal
dashed line in panel (c) indicates the ratio (N h I /N h B ) for pure water.

the first peak height and the average value of N h at the bulk
region. Clearly, the ratio increases with x G . This suggests thatthe increased number of glycerol molecules draw more watermolecules at the interfacial region from the outer solvation
layers. Hence, the ratio NhI /NhB increases with increasing glycerol concentration in the mixture. However, we have alreadywitnessed in Fig. 3 that relative “tetrahedrality” (hqI i/hqB i)of solvent remains almost insensitive to x G . Thus glycerol,on the one hand, uniformly disrupts the local tetrahedral
H-bonding structure of water at both the bulk and interface
and, on the other hand, it retains some water· · · water H-bondsin the solute’s first hydration shell. This is an excellent scenario
for stabilization of protein at low temperature.
The overall decrease of N h with increasing DMSO concentration in the DW mixture is evidenced in Fig. 6(b). However, the situation is substantially different here comparedto the GW mixture. As seen in Fig. 6(c), the ratio (NhI /NhB )FIG. 7. Density of water as a function of distance from the solute for 4 different compositions of the two binary mixtures: (a) GW and (b) DW. The results
are shown here only for 4 compositions of each of the mixtures.

134501-8

Snehasis Daschakraborty

While ρw peaks at r ∼ 4.0 Å distance from the solute, the firstminimum arises at r ∼ 5.5 Å. The density almost levels off
after 10 Å distance from the solute, which is the bulk regime.
With increasing concentration of glycerol and DMSO in the
GW and the DW mixture, respectively, the density of water
decreases at all regions varying the distance from the solute.
This is pretty straightforward since the total number of watermolecules decreases with increasing glycerol or DMSO concentration. However, the question remains whether the densityof water decreases uniformly at all regions around the solute.
We discuss this now.
Figure 8(a) exhibits the density of the interfacial water
solvent ρI (the first peak intensities) and the bulk waterρB as a function of glycerol and DMSO concentration for
GW and DW binary mixtures, respectively. For GW mixture, while ρI decreases from ∼2.2 (for neat water) to ∼1.6 (atx G = 0.15), ρB decreases from ∼1.0 (for neat water) to ∼0.6
(at x G = 0.15). Therefore, the reduction of water density at theinterface (∼27%) is less than that at the bulk (∼40%). On theother hand for the DW mixture, while the first peak intensity
decreases from ∼2.2 (for neat water) to ∼1.2 (at x D = 0.15),
the bulk water density decreases from ∼1.0 (for neat water)to ∼0.6 (at x D = 0.15). Therefore, the interfacial water solvent
density decreases by 45% for the DW mixture at x D = 0.15.
This reduction is almost twice the reduction of the interfacial
water solvent density for the same increase in the glycerol
concentration.
Figure 8(b) plots the ratio ρI /ρB as a function of theconcentration of glycerol and DMSO for the GW and DW
binary mixtures, respectively. Interestingly, the ratio increases
with increasing x G for GW binary mixtures. In contrast, theFIG. 8. (a) The densities of the interfacial water molecules (ρI ) and the bulk
water molecules (ρB ) as a function of glycerol and DMSO concentration for
GW (blue dashed lines) and DW (red solid lines) binary mixture, respectively.
(b) The ratio, ρI /ρB , as a function of the composition of the two mixtures.
The dashed lines in panel (b) are the linear fits of the data.

J. Chem. Phys. 148, 134501 (2018)

ratio decreases with DMSO concentration in the DW mixture. This provides an absolutely clear picture that glycerol
attracts solvent water molecules from the bulk region at ahigher concentration of glycerol and thereby prevents complete drying out of water from the first hydration shell of
the hydrophobic solute. In addition, glycerol acts as a goodanti-freezing substance by disrupting both the local tetrahedral structure and water· · · water H-bonding network. Glycerol, therefore, endorses two level protections of living cells
and biomolecules at the extreme condition. On the other
hand, DMSO excludes water molecules from the first hydration shell at higher concentration of DMSO. These are consistent with the changes of the H-bond numbers per watermolecule in these mixtures as a function of cosolvent concentration. Clearly, at high concentration of DMSO, water density
becomes significantly low near the solute and therefore DMSO
is not a very good cosolvent for stabilizing biomolecules atthe extreme condition particularly at a high concentrationof the DMSO mixture. This is consistent with experimental
observation.58
D. Hydrophobic effect and orientational orderingof the cosolvents

In Secs. III A–III C, we observed the contrasting roles ofglycerol and DMSO on the structural disruption of water solvent’s orientational and H-bond structure at the hydrophobic
interface and in the bulk region. While glycerol affects the regular tetrahedral H-bonding structure of water almost uniformly
at all regions, DMSO’s role is somewhat non-uniform since thelatter affects the interfacial water structure more intensely thanthe bulk water. These results corroborate with the observations
from earlier studies. In this section, we focus on the plausible
mechanism for the above contrasting behaviors of these twocosolvents. We see whether the contrasting behavior is somehow correlated with the intrinsic molecular structures of thesetwo cosolvents.
Here, we examine the orientational distribution of the key
segments of the two cosolvents at different distances fromthe hydrophobic solute. Figure 9(b) exhibits two-dimensional
probability density distribution (P(r,θ)) of the angles
θ—between glycerol CT . . . Xe vector and the glycerol CG −−CTbond vector [see Fig. 9(a) for the atom labels]—and the distance rXeCT between the Xe atom and the glycerol CT atom
for the GW mixture with x G = 0.10, as a representative composition. The distribution has a high-intensity peak at rXeCT= 3.5 Å and θ = 0◦ angle. Therefore, the terminal CH2
groups align toward the hydrophobic solute. This is due to thehydrophobic interaction between the CH2 groups of glycerol
and the Xe solute.
Now, we study the orientational alignment of the hydroxyl
groups of the glycerol molecules. We have plotted—in
Fig. 9(c)—P(r,φ) as functions of the angle φ {between theglycerol OH−−H bond vector and glycerol OH. . . Xe vector [seeFig. 9(a) for the atom labels]} and the distance rXeOH (between
the Xe atom and the OH atom of the glycerol molecule) for
the GW mixture with x G = 0.10. Note that similar pictures (notshown here) are seen for other compositions. In the interfacial
region (rXeOH < 5.5 Å), the distribution has two maxima:
one at φ = 70◦ and the other one at φ = 180◦ . Therefore,

134501-9

Snehasis Daschakraborty

J. Chem. Phys. 148, 134501 (2018)

Fig. S5 of the supplementary material.) This water-like orientational arrangement is the possible factor, which does not allow
the glycerol molecules to discriminate between the interfacial
and the bulk water molecules while perturbing the tetrahedral
H-bond structure of water.
Now we turn our focus on the orientational arrangementof DMSO molecules around the solute. Figure 10(b) exhibits
two-dimensional probability density distribution (P(r, θ)) of
the angle θ—between the DMSO CD−−Xe vector and DMSO
S−−CD bond vector—and the distance rXeCD —between the
Xe atom and DMSO CD atom—for the DW mixture withx D = 0.10, as a representative one. rXeCD and θ are labeled inFig. 10(a). Note that similar pictures (not shown here) are seen

FIG. 9. (a) Schematic diagrams of a glycerol molecule and the Xe solute.
The angle θ between the CT −−Xe vector and the CG −−CT bond vector andthe angle φ between the glycerol OH−−H vector and glycerol OH−−Xe vector
are indicated. (b) Two-dimensional probability density distribution (P(r,θ))of the angle θ and the distance between Xe and the CT atom of the glycerol
molecule (rXeC T ) for the GW mixture at a representative mole fraction x G= 0.10. (c) Two-dimensional probability density distribution (P(r,φ)) of theangle φ and the distance between Xe and the O atom of the glycerol molecule(r XeOH ) for the GW mixture at a representative mole fraction x G = 0.10. The
vertical dashed lines in all panels indicate the interfacial and bulk regions.

the majority of the OH groups of the glycerol molecules areeither aligned tangentially to the hydrophobic solute surface
or they are aligned toward the bulk region. Apparently, thisangular arrangement of glycerol’s hydroxyl groups is very
similar to the angular arrangement of water OH bonds inthe hydrophobic interfacial region. The latter is presented inFig. S4 of the supplementary material. Majority of the water
OH bonds are aligned either tangential to the solute surface
or are directed toward the bulk water molecules. Therefore,
the glycerol OH groups would participate in both tangential and bulk H-bonding with solvent water molecules. (SeeFIG. 10. (a) Schematic diagrams of a DMSO molecule and the Xe solute.
The angle θ between the CD−−Xe vector and the S−−CD bond vector and theangle φ between the DMSO S==O vector and S−−Xe vector are indicated.
(b) Two-dimensional probability density distribution (P(r,θ)) of the angle θand the distance between Xe and the CD atom of DSMO molecule (r XeCD )for the DW mixture at a representative mole fraction x D = 0.10. (c) Twodimensional probability density distribution (P(r,φ)) of the angle φ and the
distance between Xe and the O atom of DSMO molecule (r XeO ) for the DW
mixture at a representative mole fraction x D = 0.10. The vertical dashed lines
in panels (b) and (c) indicate the interfacial and bulk regions.

134501-10

Snehasis Daschakraborty

for other compositions. The peak is centered at θ = 0◦ angle
and rXeCD = 3.5 Å distance from the solute. Therefore—similar
to the alignment of terminal CH2 groups of interfacial glycerol molecules—the methyl groups of DMSO are also directed
toward the hydrophobic Xe solute. This gives a strong indication of hydrophobic interaction between the DMSO’s methylgroups and the Xe solute.
Finally, we focus on the angular alignment of the S==O
groups of the DMSO cosolvent molecules in the proximity of the solute. Figure 10(c) plots P(r,φ) as a function of
the angle φ—between the DMSO S. . . Xe vector and DMSOS==O bond vector [see Fig. 10(a) for the labels]—and the distance rXeO —between the Xe atom and DMSO O atom—for
the DW mixture having x D = 0.10, as a representative one.
Clearly, for the interfacial DMSO molecules P(r,φ) peaks at
φ ∼ 65◦ . This indicates that the S==O groups of the interfacial DMSO molecules are aligned tangentially to the solute.
Therefore, these S==O groups will involve in H-bonding onlywith water· · · water tangential H-bonds and thereby selectively
affect the tangential H-bonding structure of the water solventin the close proximity of the Xe solute. We have shown in
Fig. S5 of the supplementary material that the glycerol
cosolvent molecules affect both the tangential and thebulk H-bonding of interfacial water solvent molecules,
while DMSO molecules selectively affect the tangential
H-bonds. Note that these tangential water· · · water H-bonds
are a unique identity of the water hydrophobic interfacial
region. Therefore, the selective attacking of the tangential H-bonds by DMSO molecules perturbs the interfacial
hydration layer very differently than that in bulk. On theother hand, glycerol molecules behave very similar to thoseof solvent water molecules in terms of H-bonding withsolvent water. These emphasize why glycerol influences
the water solvent structure uniformly at all regions, but
DMSO has preferential influence on interfacial water solvent
structure.
IV. CONCLUSIONSIn this paper, we have systematically studied the effect of
two well-known cryoprotectants—glycerol and DMSO—on
the local tetrahedral order of water molecules at 5 ◦ C adjacent to a hydrophobic solute and at the bulk region. We havefound a steady reduction of tetrahedral order and the numberof water· · · water H-bonds per water molecule at both the bulk
and hydrophobic interface with increasing concentration ofglycerol and DMSO in water. Contrasting roles of glycerol and
DMSO have been clearly evidenced. While the interfacial andbulk tetrahedral H-bonding networks of water are affected byglycerol almost uniformly, DMSO disrupts the water structurearound the nonpolar solute in a nonhomogeneous way. DMSO
selectively affects the interfacial water structure more than thatof the bulk water. Also, the presence of glycerol cosolvent in
the mixture does not allow a complete drying out of the interface by retaining some water molecules in this region. Notethat this cosolvent induced water retaining is the key for protein stabilization against denaturation. In contrast, the strongerhydrophobic interaction between hydrophobic methyl residuesof DMSO and the hydrophobic solute replaces water molecules

J. Chem. Phys. 148, 134501 (2018)from the interface. This provides somewhat clear indications
that glycerol has more efficacy for stabilizing protein than
DMSO in aqueous solution. This is consistent with previous
studies.56–64 The contrasting role of glycerol and DMSO on
disrupting the local order of interfacial water has been connected with hydrophobic interaction (between the solute andthe cosolvents) induced orientation of the cosolvents adjacentto the hydrophobic solute. While glycerol has three hydroxyl
groups, forming strong HBs with solvent water by both accepting and donating H-bonds, DMSO can only accept the H-bondfrom water solvents. Therefore, the glycerol molecules behave
almost like water solvent (in terms of H-bonding structure) andtherefore can act as water replacement with minimal distortionof the overall H-bond network near the hydrophobic interface.
Finally, this work has mainly focused on the disruptionof the tetrahedral structure of liquid water at the hydrophobic interface relative to the bulk. However, it is also essential
to know how these cosolvents affect the dynamics of watermolecules in different regions from the solute. The results ofthis work will be presented separately elsewhere.
SUPPLEMENTARY MATERIAL

See supplementary material for (i) a table detailing compositions and simulated densities of different mixtures, (ii)solute-solvent radial distribution function, (iii) plot of the number of H-bonds per water molecule as a function of distancefrom the solute, (iv) two-dimensional probability density distribution (P(r,θ)) of the angles θ—between the water OW . . . Xe
vector and the water OW −−HW bond vector—and the distancerXeCT —between the Xe atom and the glycerol CT atom for the
GW mixture—with x G = 0.10, and (v) angular distribution ofcosolvent-water H-bond at the interface.
ACKNOWLEDGMENTS

I acknowledge Dr. Sandipa Indra for various suggestions.
I acknowledge the Indian Institute of Technology Patna forcomputational facilities and financial assistance.
The author declares no competing financial interest.
1 J. L. Dashnau, N. V. Nucci, K. A. Sharp, and J. M. Vanderkooi, J. Phys.

Chem. B 110, 13670 (2006).
2 J. J. Towey, A. K. Soper, and L. Dougan, J. Phys. Chem. B 120, 4439

(2016).
3 J. J. Towey, A. K. Soper, and L. Dougan, J. Phys. Chem. B 116, 13898(2012).
4 J. J. Towey and L. Dougan, J. Phys. Chem. B 116, 1633 (2012).
5 D.

B. Wong, K. P. Sokolowsky, M. I. El-Barghouthi, E. E. Fenn,
C. H. Giammanco, A. L. Sturlaugson, and M. D. Fayer, J. Phys. Chem.
B 116, 5479 (2012).
6 B. Kirchner and M. Reiher, J. Am. Chem. Soc. 124, 6206 (2002).
7 K.-I. Oh, K. Rajesh, J. F. Stanton, and C. R. Baiz, Angew. Chem., Int. Ed.
56, 11375 (2017).
8 A. K. Soper and A. Luzar, J. Chem. Phys. 97, 1320 (1992).
9 A. K. Soper and A. Luzar, J. Phys. Chem. 100, 1357 (1996).
10 M. R. Harpham, N. E. Levinger, and B. M. Ladanyi, J. Phys. Chem. B 112,283 (2008).
11 C. Chen, W. Z. Li, Y. C. Song, and J. Yang, J. Mol. Liq. 146, 23 (2009).
12 L. Weng, C. Chen, J. Zuo, and W. Li, J. Phys. Chem. A 115, 4729 (2011).
13 J. Bachler, V. Fuentes-Landete, D. A. Jahn, J. Wong, N. Giovambattista,
and T. Loerting, Phys. Chem. Chem. Phys. 18, 11058 (2016).
14 A. V. Egorov, A. P. Lyubartsev, and A. Laaksonen, J. Phys. Chem. B 115,
14572 (2011).

134501-11

Snehasis Daschakraborty

15 N. Engel, K. Atak, K. M. Lange, M. Gotz, M. Soldatov, R. Golnak,

E. Suljoti, J.-E. Rubensson, and E. F. Aziz, J. Phys. Chem. Lett. 3, 3697
(2012).
16 A. Perera and R. Mazighi, J. Chem. Phys. 143, 154502 (2015).
17 S. Indra and R. Biswas, J. Phys. Chem. B 120, 11214 (2016).
18 P. B. Conrad and J. J. de Pablo, J. Phys. Chem. A 103, 4049 (1999).
19 I. I. Vaismant and M. L. Berkowitz, J. Am. Chem. Soc. 114, 7889 (1992).
20 S. Chowdhuri and A. Chandra, J. Chem. Phys. 119, 4360 (2003).
21 A. Luzar and D. Chandler, J. Chem. Phys. 98, 8160 (1993).
22 C. Branden and J. Tooze, Introduction to Protein Structures, 2nd ed. (Taylor
and Francis, London, 1999).
23 T. E. Creighton, Proteins, 2nd ed. (Freeman, New York, 1999).
24 C. O. Fagain, Protein Stability and Stabilization of Protein Function(Landes Bioscience, Georgetown, TX, 1997).
25 S. A. Hawley, Biochemistry 10, 2436 (1971).
26 P. L. Privalov, J. Chem. Thermodyn. 29, 447 (1997).
27 R. Kitahara, A. Okuno, M. Kato, Y. Taniguchi, S. Yokoyama, andK. Akasaka, Magn. Reson. Chem. 44, S108–S113 (2006).
28 D. Sanfelice and P. A. Temussi, Biophys. J. 208, 4 (2016).
29 P. L. Privalov, Crit. Rev. Biochem. Mol. Biol. 25, 281 (1990).
30 A. Pastore, S. R. Martin, A. Politou, K. C. Kondapalli, T. Stemmler, and
P. A. Temussi, J. Am. Chem. Soc. 129, 5374 (2007).
31 V. R. Agashe and J. B. Udgaonkar, Biochemistry 34, 3286 (1995).
32 P. L. Privalov, Y. V. Griko, S. Y. Venyaminov, and V. P. Kutyshenko, J.
Mol. Biol. 190, 487 (1986).
33 C. Yang, S. Jang, and Y. Pak, Nat. Commun. 5, 5773 (2014).
34 C. L. Dias, T. Ala-Nissila, J. Wong-ekkabut, I. Vattulainen, M. Grant, and
M. Karttunen, Cryobiology 60, 91 (2010).
35 F. Franks, R. H. M. Hatley, and H. L. Friedman, Biophys. Chem. 31, 307(1988).
36 F. Franks, Adv. Protein Chem. 46, 105 (1995).
37 B. L. Chen and J. A. Schellman, Biochemistry 28, 685 (1989).
38 B. Hallerbach and H. J. Hinz, Biophys. Chem. 76, 219 (1999).
39 C. L. Dias, T. Ala-Nissila, M. Karttunen, I. Vattulainen, and M. Grant,
Phys. Rev. Lett. 100, 118101 (2008).
40 M. Marqués, J. M. Borreguero, H. E. Stanley, and N. V. Dokholyan, Phys.
Rev. Lett. 91, 138103 (2003).
41 M. Davidovic, C. Mattea, J. Qvist, and B. Halle, J. Am. Chem. Soc. 131,
1025 (2009).
42 D. Sanfelice, E. Morandi, A. Pastore, N. Niccolai, and P. A. Temussi,
ChemPhysChem 16, 3599 (2015).
43 C. J. Tsai, J. V. Maizel, Jr., and R. Nussinov, Crit. Rev. Biochem. Mol.
Biol. 37, 55 (2002).
44 C. F. Lopez, R. K. Darst, and P. J. Rossky, J. Phys. Chem. B 112, 5961
(2008).
45 J. R. Grigera and A. N. McCarthy, Biophys. J. 98, 1626 (2010).
46 E. Ascolese and G. Graziano, Chem. Phys. Lett. 467, 150 (2008).
47 M. Maiti, S. Weiner, S. V. Buldyrev, H. E. Stanley, and S. Sastry, J. Chem.
Phys. 136, 044512 (2012).
48 C. L. Dias, Phys. Rev. Lett. 109, 048104 (2012).
49 S. Parui and B. Jana, J. Phys. Chem. B 121, 7016 (2017).
50 S. B. Kim, J. C. Palmer, and P. G. Debendetti, Proc. Natl. Acad. Sci. U. S. A.
113, 8991 (2016).
51 T. Hajari and S. Bandyopadhyay, J. Chem. Phys. 146, 225104 (2017).
52 A. Ben-Naim, Adv. Biol. Chem. 3, 29 (2013).
53 V. Bianco and G. Franzese, Phys. Rev. Lett. 115, 108101 (2015).
54 M. E. Johnson, C. Malardier-Jugroot, and T. Head-Gordon, Phys. Chem.
Chem. Phys. 12, 393 (2010).
55 A. M. Schrader, C.-Y. Cheng, J. N. Israelachvili, and S. Han, J. Chem.
Phys. 145, 041101 (2016).
56 K. Gekko and S. N. Timasheff, Biochemistry 20, 4667 (1981).
57 V. Vagenende, M. G. S. Yap, and B. L. Trout, Biochemistry 48, 11084(2009).
58 S. Bhattacharjya and P. Balaram, Proteins 29, 492 (1997).
59 K. Hamaguchi, J. Biochem. 56, 441 (1964).
60 P. Huang, A. Dong, and W. S. Caughey, J. Pharmacol. Sci. 84, 387 (1995).
61 T. Kamiyama, T. Matsusita, and T. Kimura, J. Chem. Eng. Data 48, 1301
(2003).
62 S. Roy, B. Jana, and B. Bagchi, J. Chem. Phys. 136, 115103 (2012).
63 S. Roy and B. Bagchi, J. Phys. Chem. B 117, 4488 (2012).
64 S. Roy and B. Bagchi, J. Chem. Phys. 139, 034308 (2013).
65 D. T. Bowron, A. Filipponi, M. A. Roberts, and J. L. Finney, Phys. Rev.
Lett. 81, 4164 (1998).

J. Chem. Phys. 148, 134501 (2018)
66 H. S. Ashbaugh, D. Asthagiri, L. R. Pratt, and S. B. Rempe, Biophys.

Chem. 105, 323 (2003).
67 N. Galamba, J. Phys. Chem. B 117, 2153 (2013).
68 N. Galamba, J. Phys. Chem. B 118, 4169 (2014).
69 N. Galamba, J. Phys. Chem. B 117, 589 (2013).
70 B. Lee, Biopolymers 24, 813 (1985).
71 R. L. Baldwin, Proc. Natl. Acad. Sci. U. S. A. 111, 13052 (2014).
72 J. T. Titantah and M. Karttunen, J. Am. Chem. Soc. 134, 9362 (2012).
73 J. Grdadolnika, F. Merzela, and F. Avbelj, Proc. Natl. Acad. Sci. U. S. A.

114, 322 (2017).
74 V. I. Artyukhov, A. Yu. Pulver, A. Peregudov, and I. Artyuhov, J. Chem.

Phys. 141, 034503 (2014).
75 R. Haselmeier, M. Holz, W. Marbach, and H. Weingärtner, J. Phys. Chem.

99, 2243 (1995).
76 D. Paschek, J. Phys. Chem. 120, 6674 (2004).
77 O. Coskuner and U. K. Dieters, Z. Phys. Chem. 221, 785 (2007).
78 K. Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 103, 4570 (1999).
79 D. M. Huang and D. Chandler, Proc. Natl. Acad. Sci. U. S. A. 97, 8324

(2000).
80 H. S. Ashbaugh and M. E. Paulaitis, J. Am. Chem. Soc. 123, 10721(2001).
81 S. Rajamani, T. M. Truskett, and S. Garde, Proc. Natl. Acad. Sci. U. S. A.

102, 9475 (2005).
82 H. S. Ashbaugh and L. R. Pratt, Rev. Mod. Phys. 78, 159 (2006).
83 M. V. Athawale, S. N. Jamadagni, and S. Garde, J. Chem. Phys. 131, 115102(2009).
84 S. Garde and A. J. Patel, Proc. Natl. Acad. Sci. U. S. A. 108, 16491(2011).
85 C. Camilloni, D. Bonetti, A. Morrone, R. Giri, C. M. Dobson, M. Brunori,

S. Gianni, and M. Vendruscolo, Sci. Rep. 6, 28285 (2016).
86 R. V. Rariy and A. M. Klibanov, Proc. Natl. Acad. Sci. U. S. A. 94, 13520(1997).
87 T. Arakawa, Y. Kita, and S. N. Timasheff, Biophys. Chem. 131, 62(2007).
88 J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J.

Comput. Chem. 25, 1157 (2004).
89 R. Chelli, P. Procacci, G. Cardini, R. G. DellaValle, and S. Califano, Phys.

Chem. Chem. Phys. 1, 871 (1999).
90 J. Blieck, F. Affouard, P. Bordat, A. Lerbret, and M. Descamps, Chem.

Phys. 317, 253 (2005).
91 J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).
92 F. O. Akinkunmi, D. A. Jahn, and N. Giovambattista, J. Phys. Chem. B

119, 6250 (2015).
93 M. L. Strader and S. E. Feller, J. Phys. Chem. A 106, 1074 (2002).
94 A. D. Becke, Phys. Rev. A 38, 3098 (1988).
95 C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
96 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb,

J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, and G. A.
Petersson, gaussian 09, Revision B.01, Gaussian, Inc., Wallingford, CT,2010.
97 L. Martı́nez, R. E. Andrade, G. Birgin, and J. M. Martı́nez, J. Comput.
Chem. 30, 2157 (2009).
98 D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and
H. J. C. Berendsen, J. Comput. Chem. 26, 1701 (2005).
99 S. A. Nosé, J. Chem. Phys. 81, 511 (1984).
100 W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
101 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, andJ. R. Haak, J. Chem. Phys. 81, 3684 (1984).
102 M. R. Shirts, J. Chem. Theory Comput. 9, 909 (2013).
103 M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
104 S. A. Schichman and R. A. Amey, J. Phys. Chem. 75, 98 (1971).
105 N. Giovambattista, P. G. Debenedetti, F. Sciortino, and H. E. Stanley, Phys.
Rev. E 71, 061505 (2005).
106 Z. Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. Debenedetti, and
H. Stanley, Phys. Rev. E 76, 051201 (2007).
107 S. Chatterjee, P. G. Debenedetti, F. H. Stillinger, and R. M. Lynden-Bell,J. Chem. Phys. 128, 124511 (2008).
108 D. R. Nutt and J. C. Smith, J. Am. Chem. Soc. 130, 13066 (2008).
109 M. Paolantoni, N. F. Lago, M. Albertı́, and M. A. Lagana, J. Phys. Chem.
A 113, 15100 (2009).
110 E. Dubouev́-Dijon and D. Laage, J. Phys. Chem. B 119, 8406
(2015).
111 S. Indra and S. Daschakraborty, Chem. Phys. Lett. 685, 322 (2017).