Entropic enrichment of cosolvent near a very large solute immersed in solvent-cosolvent binary mixture: Anomalous dependence on bulk cosolvent concentration

Journal of Molecular Liquids 247 (2017) 403–410

Contents lists available at ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

Entropic enrichment of cosolvent near a very large solute immersed in solvent-cosolvent binary mixture: Anomalous dependence on bulk cosolvent concentration Masahiro Kinoshita ⁎, Tomohiko Hayashi Institute of Advanced Energy, Kyoto University, Uji, Kyoto 611-0011, Japan

a r t i c l e

i n f o

Article history: Received 7 August 2017 Received in revised form 25 September 2017 Accepted 26 September 2017 Available online 28 September 2017 Keywords: Cosolvent Excluded volume Entropic potential Enrichment Integral equation theory

a b s t r a c t We analyze the number density profiles of solvent and cosolvent entropically formed near a very large solute using an integral equation theory combined with a reliable bridge function. The entropic enrichment of cosolvent near the solute exhibits rather complex behavior, depending on the solvent and cosolvent sizes and the cosolvent concentration in the bulk. For a system where the cosolvent particles are sufficiently larger than the solvent particles, when the cosolvent concentration in the bulk XC is changed as a parameter, the number of cosolvent particles in contact with the solute NContact becomes largest at a surprisingly low value of XC: XC ~ 0.01 for dC/dS = 4, ηT = 0.3831, and dU/dS = 20 (dU, dS, and dC are the particle diameters of the solute, solvent, and cosolvent, respectively, and ηT is the total packing fraction of the solvent-cosolvent mixture). The largest value of NContact is approximately twice larger than NContact in pure cosolvent. The bulk solvent-cosolvent mixture is thermodynamically stable as a single phase, and the intriguing behavior mentioned above is induced by the introduction of solute surface into the solvent-cosolvent mixture. © 2017 Elsevier B.V. All rights reserved.

1. Introduction In general, the structure and properties of fluid near a solid surface or a large solute (i.e., surface-induced structure and properties) often differ from those of the bulk fluid. In particular, when the fluid comprises multiple components, the fluid composition within the surface-induced layer can be remarkably different from that in the bulk. Conversely, the formation of surface-induced layer can have large effects on the polyatomic structure of the solute. An example is the destabilization or stabilization of the native structure of a protein caused by adding a cosolvent to water [1–8]. The investigation of this cosolvent effect is important in the respect that it provides a clue to the folding/unfolding mechanism of a protein. Urea [1,2] and its alkylated derivatives [3] are known as typical denaturants. Whether a cosolvent is enriched near the protein surface or not is an important factor [7,8]. When it is not enriched, the structural stability can indirectly be influenced by the cosolvent through the modification of water structure in the bulk. When it is enriched, on the other hand, significantly many solvent molecules are replaced by cosolvent molecules near the protein surface, and the differences between protein-cosolvent and protein-solvent interaction energies can change the protein structure. For instance, urea (or its alkylated derivative) is significantly enriched and many urea molecules ⁎ Corresponding author. E-mail address: [email protected] (M. Kinoshita).

https://doi.org/10.1016/j.molliq.2017.09.108 0167-7322/© 2017 Elsevier B.V. All rights reserved.

come in contact with the protein surface. This enrichment lowers the structural stability for the following reasons [6,8]: the van der Waals (vdW) interaction energy for protein-urea is substantially lower than that for protein-water; and in terms of the vdW interaction energy, a structure with larger area of the protein surface exposed to the solvent-cosolvent mixture (i.e., denatured state) is more favored. A question then arises: why is urea enriched near the protein surface? It is enriched because the protein-urea affinity is higher than the waterurea one. This answer is, however, quite ambiguous since the affinities are determined by the competition of energetic and entropic factors. The principal factor causing the urea enrichment is not necessarily the protein-urea vdW interaction. We note that the entropic factor is omnipresent in any solution system but more difficult to elucidate than the energetic one when the solvent comprises multiple components (one of the components and the others are referred to as “solvent” and “cosolvents”, respectively). The degree of cosolvent enrichment near the solute is largely influenced by an entropic factor explained in what follows. The presence of a solute generates an excluded space (ES) which the centers of solvent and cosolvent molecules cannot enter. The volume of the ES is referred to as “excluded volume (EV)”. The presence of a solvent or cosolvent molecule also generates an ES for the other solvent and cosolvent molecules. When a solvent or cosolvent molecule contacts the solute, the overlap of two ESs occurs. Consequently, the total volume available to the other solvent and cosolvent molecules increases by the volume of

404

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

overlapped space, leading to a gain of their entropy. Thus, entropic solute-solvent and solute-cosolvent affinities come into play, driving solvent and cosolvent molecules to contact the solute (i.e., promoting their enrichment near the solute). This can be referred to as “entropic EV effect” [9,10]. On the other hand, it is entropically more favorable for all the solvent and cosolvent molecules to move to and fro within the whole system. The degree of cosolvent and solvent enrichment is determined by the competition of these two factors. For the investigation of entropic behavior of a solution system, the molecules constituting the system are all modeled as neutral hard bodies interacting through simple hard-body potentials with no attractive tails. In such a system, all of the allowed system configurations share the same energy and the system behavior is purely entropic in its origin. The essential physical quantities are the EV generated by a molecule for the other molecules and the total packing fraction of the system, and a sphere is often chosen for the molecular shape as the most fundamental geometry for investigating the effects of these physical quantities. The properties of bulk binary mixtures of hard spheres with different sizes have been studied rather extensively [11,12]. Also, the interaction induced between large hard spheres immersed in small hard spheres [13,14] or in binary [15] and multi-component [16] mixtures of hard spheres with different sizes has also been studied and a significant amount of information is available. On the other hand, the surface-induced structure entropically formed remains rather elusive. It is well known that cosolvent and solvent molecules are entropically enriched near a large solute [14,15]. It is apparent that the entropic EV effect becomes larger and the enrichment is enhanced as the solute size increases. However, the effect of complex interplay of the solvent and cosolvent sizes and the cosolvent concentration in the bulk solventcosolvent mixture has not been comprehended yet. In this study, we analyze the entropic enrichment of cosolvent near a very large solute immersed in a solvent-cosolvent binary mixture at infinite dilution using an integral equation theory [17], a statistical-mechanical theory for fluids. The solute, cosolvent, and solvent molecules are all modeled as neutral hard spheres. The use of this system allows us to exclusively investigate the entropic EV effect. The number densities of solvent and cosolvent at the solute surface are calculated by changing XC (XC is the cosolvent mole fraction in the bulk) from 0 to 1 with a small increment for two representative, different values of dC/ dS (dC and dS are diameters of the cosolvent and solvent particles, respectively). A counterintuitive finding is as follows: for sufficiently large dC/dS, the cosolvent number density at the solute surface becomes highest when XC (XC is the cosolvent mole fraction in the bulk) is much lower than 1: XC ~ 0.01 for dC/dS = 4, ηT = 0.3831, and dU/dS = 20 (dU is the solute diameter and ηT is the total packing fraction of the solventcosolvent mixture). It is known that a binary mixture of hard spheres with different diameters undergoes phase separation if ηT is sufficiently high and the diameter ratio exceeds a threshold value. We verify that no such phase separation occurs under the conditions considered in this study. The anomalous behavior originates not from the solvent-solvent, solvent-cosolvent, and cosolvent-cosolvent correlations but from the interplay of not only these three correlations but also the solute-solvent and solute-cosolvent correlations.

2. Model and theory 2.1. Solute, solvent, and cosolvent models The solute is immersed in a solvent-cosolvent binary mixture at infinite dilution. The solvent and cosolvent particles are modeled as neutral hard spheres with diameters dS and dC (dS b dC), respectively (the subscript “S” and “C” denote “solvent” and “cosolvent”, respectively). The solute is also modeled as a neutral hard sphere whose diameter is dU. All of the hard spheres interact through hard-sphere potentials with no attractive tails.

We consider system 1 (dC = 2dS) and system 2 (dC = 4dS) where dU is set at 20dS. That is, we consider two cosolvents with dC = 2dS and 4dS, respectively. The solvent and cosolvent number densities in the bulk are denoted by ρS and ρC, respectively, and ρCd3C is changed from 0.0 to 0.7317 as an important parameter under the condition that “ρSd3S + ρCd3C” is kept at 0.7317, i.e., ηT = π(ρSd3S + ρCd3C) / 6 is kept at 0.3831. The value, 0.3831, is pertinent to water at 298 K and 1 atm. What matters is to set ηT at a sufficiently high value so that the entropic effect can become substantial. Changing 0.3831 to, for example, 0.37 or 0.40 is acceptable. In addition, changing ηT as a function of the cosolvent concentration might be more realistic. However, our conclusions are not likely to be altered by these treatments. To extract the purely entropic effect from the aqueous solution system, setting ηT at the water value seems to be the most reasonable. XC = ρC / (ρS + ρC), which represents the cosolvent mole fraction in the bulk, is referred to as “cosolvent concentration”. “XS = 1 − XC” is the solvent mole fraction in the bulk. XC is changed from 0.0 to 1.0 with a small increment. We note that in this study the component with a larger particle diameter is categorized as the cosolvent and cases of XS b XC are also analyzed to investigate the basic behavior of the entropic EV effect. 2.2. Integral equation theory We employ an integral equation theory (IET) for simple fluids [17]. It comprises the Ornstein-Zernike (OZ) equation and a closure equation. The calculation procedure consists of the following two steps: (1) the solvent-solvent, solvent-cosolvent, and cosolvent-cosolvent correlation functions are calculated; (2) and the solute-solvent and solutecosolvent correlation functions are calculated using the correlation functions obtained in step (1) as part of the input data. The hypernetted-chain (HNC) closure equation neglecting the bridge function is known to be quantitatively unreliable when XC is finite and dC/dS takes a very large value. In addition, the solute-induced density structure calculated using the HNC is rather inaccurate for a very large solute. We note that the system treated in this study features very high size asymmetry: dC/dS = 4 and dU/dS = 20. Attard and Patey [18] calculated all of the bridge diagrams with two and three field points using Monte Carlo integration techniques and included those with more than three field points by the Páde approximant. It could be applicable to the present system in principle but very time-consuming. Several closure equations other than the HNC or different types of bridge functions which can be calculated through much more practical routes were proposed [19–24]. They were tested for bulk pure fluids and bulk mixtures with low size asymmetry and shown to be sufficiently accurate. However, it is not definite if they can successfully be applied to a system with high size asymmetry like the present one. In fact, for some of the closure equations and bridge functions, we showed that they become problematic and give pathological results (i.e., results which are much worse than those from the HNC) [25]. (Also, it is well known that the PercusYevick (PY) closure is better than the HNC for a pure hard-sphere solvent [17] but much worse when a large solute is immersed in the solvent [18] and that the PY is inferior to the HNC for a binary mixture of hard spheres unless its size asymmetry is low.) Hence, we adopt the functional form for the bridge function proposed by Kinoshita [25]. In Kinoshita's functional form, the bridge function bij(r) is expressed as n o2 n o bij ðr Þ ¼ −0:5 γ ij ðr Þ = 1 þ 0:8γ ij ðr Þ for γij ðr ÞN0;

ð1aÞ

n o2 n o bij ðr Þ ¼ −0:5 γ ij ðr Þ = 1−0:8γij ðr Þ for γ ij ðr Þb0:

ð1bÞ

Here, r is the distance between centers of two particles, γij(r) = hij(r) − cij(r), h and c, respectively, are the total and direct correlation functions, and the subscript “ij” denote the “solvent-solvent”, “solventcosolvent”, or “cosolvent-cosolvent” pair in step (1) and the “solute-

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

solvent” or “solute-cosolvent” pair in step (2). Eqs. (1a) and (1b) and its first and second derivatives with respect to γij are continuous at γij = 0. As shown in Section 3.5, the IET using Eqs. (1a) and (1b) satisfies the thermodynamic consistency [17] almost perfectly for the bulk solventcosolvent mixture. In general, when an IET is thermodynamically consistent, the result from the theory is in quantitatively good accord with that from a Monte Carlo or molecular dynamics simulation. The high reliability of Eqs. (1a) and (1b) has been demonstrated for systems with high size asymmetry and quite successful in solving several important problems including those for very large solutes immersed in multicomponent solvents [11,16,25]. The HNC possesses the following two advantages: first, it satisfies the virial-energy consistency as shown by Giacometti and coworkers [26]; second, the analytical expression for the excess chemical potential of a solute, the Morita-Hiroike formula [27,28], can be derived, which avoids the use of the direct thermodynamic integration. When the particles interact through hard-sphere potentials, the thermodynamic quantities are independent of the temperature, with the result that the virial-energy consistency is satisfied irrespective of the closure equation employed. When Eqs. (1a) and (1b) are incorporated in the closure equation, the second advantage is lost. In this study, however, the important quantities which must be calculated with sufficient accuracy are the density profiles of solvent and cosolvent near a very large solute (i.e., the surface-induced density structure of a solvent-cosolvent mixture). Therefore, the incorporation of Eqs. (1a) and (1b) is absolutely required. For numerical solution of the basic equations, a sufficiently long range R is divided into N grid points (rn = nδr, n = 0, 1, …, N − 1; δr = R / N) and all the functions are represented by their values on these points. N and δr are set at 4096 and 0.01dS, respectively. The large set of nonlinear simultaneous equations (the number of equations is 8192 in step (2)) is solved using the robust and highly efficient algorithm developed by Kinoshita and coworkers [29,30]. 3. Results and discussion The key functions which we look at are the solvent and cosolvent reduced number density profiles near the solute, which are denoted by gS(r) and gC(r), respectively. The value of gS(r) or gC(r) at contact, GS = gS(dUS) (dUS = (dU + dS) / 2) or GC = gC(dUC) (dUC = (dU + dC) / 2), represents the degree of solvent or cosolvent enrichment near the solute.

405

enriched at the solute surface. However, the degree of enrichment of the cosolvent is much higher than that of the solvent. In particular, the cosolvent enrichment is remarkable in system 2. It is interesting to note that the cosolvent number density at the solute surface ρCd3SGC in system 2 becomes highest at ρCd3S ~ 0.05 which corresponds to XC ~ 0.01. We first discuss the mixture where the cosolvent or solvent is present at infinite dilution. There are four mixtures: ρSd3S = 0.7317 and ρCd3S = 0.0 in system 1 (mixture S1); ρSd3S = 0.0 and ρCd3S = 0.09146 (ρCd3C = 0.7317) in system 1 (mixture C1); ρSd3S = 0.7317 and ρCd3S = 0.0 in system 2 (mixture S2); and ρSd3S = 0.0 and ρCd3S = 0.01143 (ρCd3C = 0.7317) in system 2 (mixture C2). Fig. 1 shows gS(r) and gC(r) in the four mixtures. GS and GC in the four mixtures are given in Table 3. We note that gS(r) and GS in mixture S1 or S2 are identical to those in pure solvent and that gC(r) and GC in mixture C1 or C2 are identical to those in pure cosolvent. When one of the cosolvents is added to the solvent at very low concentration (at infinite dilution), the cosolvent is largely enriched near the solute: the degree of enrichment is much higher than in pure cosolvent, and this trend becomes stronger for larger dC. By contrast, when the solvent is added to one of the cosolvents at very low concentration, the solvent is enriched near the solute only moderately: the degree of enrichment is lower than in pure solvent, and this trend becomes stronger for larger dC. These results can be understood by applying the Asakura-Oosawa (AO) theory [31,32] as discussed in Section 3.2. 3.2. Application of Asakura-Oosawa theory We consider a solute particle with diameter dU and another solute particle with diameter d2 immersed in particles of component 1 with diameter d1 (see Fig. 2). When component 1 is the solvent, “another solute particle” corresponds to a cosolvent particle (mixture S1 or S2). When component 1 is one the cosolvents, “another solute particle” corresponds to a solvent particle (mixture C1 or C2). The presence of solute particles generates ESs which the centers of particles of component 1 cannot enter. When the solute particles contact each other, the two ESs overlap, and the total EV decreases by the volume of overlapped space Δ Vex. By simple but tedious mathematical manipulations, we can show that ΔVex is given by 3

ΔV ex =d1 ¼ ðπ=6Þ½1 þ 3d2 dU =fd1 ðdU þ d2 Þg:

ð2Þ

3.1. Degree of enrichment of solvent or cosolvent near the solute GS, GC, ρSd3S GS, and ρCd3S GC calculated with the prescribed values of ρSd3S and ρCd3S are presented for system 1 in Table 1 and for system 2 in Table 2. It is observed that the solvent and the cosolvent are always Table 1 GS, GC, ρSGS, and ρCGC calculated with the prescribed values of ρSd3S and ρCd3S (“ρSd3S + ρCd3C” is kept at 0.7317) for system 1. ρSd3S

ρCd3S

GS

GC

ρSd3S GS

ρCd3S GC

0.7317 0.6917 0.6517 0.5717 0.4917 0.4117 0.3317 0.2517 0.1717 0.0917 0.0517 0.0000

0.00000 0.00500 0.01000 0.02000 0.03000 0.04000 0.05000 0.06000 0.07000 0.08000 0.08500 0.09146

5.719 5.558 5.402 5.099 4.809 4.533 4.268 4.015 3.774 3.544 3.433 3.295

14.02 13.37 12.75 11.57 10.49 9.488 8.567 7.718 6.937 6.220 5.885 5.475

4.184 3.845 3.520 2.915 2.365 1.866 1.416 1.011 0.648 0.325 0.178 0.000

0.0000 0.0669 0.1275 0.2314 0.3146 0.3795 0.4283 0.4631 0.4856 0.4976 0.5002 0.5008

Here, dS, dC, and dU are solvent, cosolvent, and solute diameters, respectively, and ρS and ρC are solvent and cosolvent number densities in the bulk, respectively. GS = gS(dUS) and GC = gC(dUC) where gS(r) and gC(r) denote the solvent and cosolvent reduced number density profiles near the solute, respectively.

Table 2 GS, GC, ρSGS, and ρCGC calculated with the prescribed values of ρSd3S and ρCd3S (“ρSd3S + ρCd3C” is kept at 0.7317) for system 2. ρSd3S

ρCd3S

GS

GC

ρSd3S GS

ρCd3S GC

0.73170 0.69970 0.66770 0.60370 0.53970 0.47570 0.41170 0.34770 0.28370 0.21970 0.15570 0.12370 0.09170 0.05970 0.04050 0.02898 0.01618 0.00914 0.00000

0.00000 0.00050 0.00100 0.00200 0.00300 0.00400 0.00500 0.00600 0.00700 0.00800 0.00900 0.00950 0.01000 0.01050 0.01080 0.01098 0.01118 0.01129 0.01143

5.719 5.519 5.328 4.964 4.624 4.304 4.003 3.718 3.449 3.192 2.947 2.828 2.712 2.599 2.532 2.493 2.450 2.426 2.395

60.77 55.02 49.86 41.00 33.68 27.56 22.45 18.17 14.62 11.68 9.268 8.227 7.285 6.432 5.961 5.692 5.404 5.252 5.058

4.184 3.862 3.557 2.997 2.496 2.048 1.648 1.293 0.978 0.701 0.459 0.350 0.249 0.155 0.103 0.072 0.040 0.022 0.000

0.0000 0.0275 0.0499 0.0820 0.1010 0.1103 0.1122 0.1090 0.1023 0.0935 0.0834 0.0782 0.0729 0.0675 0.0644 0.0625 0.0604 0.0593 0.0578

See Table 1 for the notations.

406

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

Fig. 1. Reduced number density profiles of solvent and cosolvent near the solute, gS(r) and gC(r), respectively. (a) gS(r) in system 1 with ρSd3S = 0.7317 and ρCd3S = 0.0 (mixture S1); gS(r) in system 2 with ρSd3S = 0.7317 and ρCd3S = 0.0 (mixture S2); these two curves are identical. (b) gC(r) in system 1 with ρSd3S = 0.7317 and ρCd3S = 0.0 (mixture S1); gC(r) in system 2 with ρSd3S = 0.7317 and ρCd3S = 0.0 (mixture S2). (c) gS(r) in system 1 with ρSd3S = 0.0 and ρCd3S = 0.09146 (mixture C1); gS(r) in system 2 with ρSd3S = 0.0 and ρCd3S = 0.01143 (mixture C2). (d) gC(r) in system 1 with ρSd3S = 0.0 and ρCd3S = 0.09146 (mixture C1); gC(r) in system 2 with ρSd3S = 0.0 and ρCd3S = 0.01143 (mixture C2). See Table 1 for the notations. The value of gC(0) for the red curve in (b) is 60.77. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In the case of d2 = dU, Eq. (2) reduces to the well known formula [13]: 3

ΔV ex =d1 ¼ ðπ=6Þf1 þ ð3=2ÞðdU =d1 Þg:

ð3Þ

Applying Eq. (2) to mixture S1 or S2 yields d1 ¼ dS ; d2 ¼ dC : :

ð4Þ

3

ΔV ex =dS ¼ ðπ=6Þ½1 þ 3dC dU =fdS ðdU þ dC Þg:

AO theory deviates significantly from the exact curve when the small spheres are dense. However, the AO theory gives a fairly accurate value of the EP at contact [13] due to fortuitous cancellation of errors (see Section 3.3): Strangely, the AO value is closer to the exact value than the value obtained by the HNC [12,14]. In these studies, only the large spheres with equal sizes are considered. Below we discuss the EP at contact estimated by the AO theory for solute spheres with different sizes for the first time. According to the AO theory [31,32], the gain of solvent entropy ΔSS for mixture S1 or S2 (ρSd3S = 0.7317) and that of cosolvent entropy

Applying Eq. (2) to mixture C1 or C2 yields d1 ¼ dC ; d2 ¼ dS : :

ð5Þ

3

ΔV ex =dC ¼ ðπ=6Þ½1 þ 3dS dU =fdC ðdU þ dS Þg:

The accuracy of the AO theory has been discussed for the whole curve of the entropic potential (EP) or force between large spheres immersed in small spheres [12,14]. It is well known that the curve from the Table 3 GS and GC calculated by the integral equation theory using the bridge function expressed by Eqs. (1a) and (1b) and comparison between (ΔSS/kB)AO and ΔSS/kB or (ΔSC/kB)AO and ΔSC/kB. Mixture

GS

GC

(ΔSS/kB)AO

(ΔSC/kB)AO

ΔSS/kB

ΔSC/kB

Error (%)

S1 C1 S2 C2

5.719 3.295 5.719 2.395

14.02 5.475 60.77 5.058

2.47 – 4.21 –

– 0.93 – 0.66

2.64 – 4.11 –

– 1.19 – 0.87

6.4 22 −2.4 24

The error is defined as Error (%) = 100{(ΔSS/kB)AO − ΔSS/kB} / (ΔSS/kB) or 100{(ΔSC/kB)AO − ΔSC/kB} / (ΔSC/kB). Mixture S1: ρSd3S = 0.7317 and ρCd3S = 0.0 in system 1. Mixture C1: ρSd3S = 0.0 and ρCd3S = 0.09146 in system 1. Mixture S2: ρSd3S = 0.7317 and ρCd3S = 0.0 in system 2. Mixture C2: ρSd3S = 0.0 and ρCd3S = 0.01143 in system 2. See Table 1 for the notations. The values of dU/d2 in mixtures S1, C1, S2, and C2 are 10, 20, 5, and 20, respectively.

Fig. 2. A solute particle with diameter dU and another solute particle with diameter d2 immersed in particles of component 1 with diameter d1. When component 1 is the solvent, “another solute particle” corresponds to a cosolvent particle. When component 1 is one of the cosolvents, “another solute particle” corresponds to a solvent particle.

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

ΔSC for mixture C1 or C2 (ρCd3C = 0.7317) are expressed as ΔSS =kB ¼ ρS dS

3

ΔSC =kB ¼ ρC dC

  3 ΔV ex =dS ;

ð6aÞ

  3 ΔV ex =dC :

ð6bÞ

3

Here, kB is Boltzmann's constant. ΔSS/kB and ΔSC/kB estimated from Eqs. (6a) and (6b) are denoted by (ΔSS/kB)AO and (ΔSC/kB)AO, respectively, and given in Table 3. ΔVex/d3S is considerably larger than ΔVex/ d3C, leading to the results mentioned at the end of Section 3.1: they can qualitatively be interpreted by the AO theory. In a quantitative sense, however, the result from the AO theory is not accurate unlike in the case of solute spheres with equal sizes, which is explained below. We note that ΔSS or ΔSC is the EP at contact: ΔSC, for example, represents the value which the EP takes when the cosolvent particle comes in contact with the solute particle. Therefore, when an exact theory is employed, the following equations hold: ΔSS =kB ¼ ln ðGC Þ;

ð7aÞ

ΔSC =kB ¼ ln ðGS Þ:

ð7bÞ

We obtain ΔSS/kB and ΔSC/kB by substituting GC and GS calculated by the IET with the reliable bridge function into Eqs. (7a) and (7b). They are compared to (Δ SS/kB)AO and (Δ SC/kB)AO, respectively, in Table 3. As shown in the table, the error is nontrivial when the AO theory is employed. Moreover, it becomes larger as dU/d2 increases. As described above, the EP at contact calculated by the AO theory is fairly accurate for spherical solutes sharing the same diameter (or for a planar wall and a spherical solute) [14]. To verify this, we consider another mixture S1 specified as ρSd3S = 0.7317, ρCd3S = 0.0, and dU = dC = 5dS. Δ SS/kB is obtained from Eq. (7a) using GC calculated by the IET with the reliable bridge function: the result is 3.325. On the other hand, (Δ SS/kB)AO estimated from Eq. (4) (dU/d2 = 1) and Eq. (6a) is 3.256. The error defined in Table 3 is only 2.1%. The AO theory certainly gives a fairly accurate estimation even in a quantitative sense when the two solute particles share the same diameter.

407

cancellation of errors and the net force calculated does not deviate seriously from the exact one. Moreover, at h = dS the AO force is zero while the exact one is repulsive, and at h = 0 the former is attractive but weaker than the latter. As a result, in the case of d2 = dU, the EP at contact, which is obtained by integrating the net force over the surface separation, can be predicted by the AO theory with sufficiently high accuracy. This success is fortuitous and no logical explanation is possible [33,34]. When dU and d2 are made different, the only change occurring is that the distribution of ρ1,Contact on the surface of one of the solute particles becomes different from that of the other. Due to this asymmetry, the fortuitous success is no more achievable and the result becomes worse as dU/d2 departs more from 1 (see Table 3). As shown by Kinoshita [33], there is only attractive component when the key and the lock modeled as a sphere and a hemispherical cavity, respectively, come very close to each other, with the result that the EP at contact estimated by the AO theory is quite erroneous. 3.4. Number density of solvent or cosolvent at the solute surface Fig. 3 shows GS and GC plotted against the bulk cosolvent concentration XC. The degree of cosolvent enrichment is significantly high in the region of low cosolvent concentration but it rapidly decreases as XC becomes higher. When the bulk mixture consists of smaller and larger particles and ηT is kept constant, the entropic EV effect becomes smaller as the concentration of larger particles increases. We define ΓS and ΓC as   3 3 Γ S ¼ ρS dS GS = ρS dS GS

Pure

;

ð8aÞ

3.3. On the physical origin of inaccuracy of Asakura-Oosawa theory We discuss the entropic force between two solute particles with diameters dU and d2 separated by h (h denotes the distance between the nearest solute surfaces), which are immersed in particles of component 1 with diameter d1. The force is induced by collisions of particles of component 1 with the solute surfaces. Let ρ1,Contact be the number density of particles of component 1 at a solute surface. For a single solute ρ1,Contact is uniform on the solute surface, but for a pair of solutes it is not uniform. The force acting at a point on the surface is in proportion to ρ1,Contact at this point. The force acts toward the center of each solute particle. As described in “Introduction”, ρ1,Contact is considerably in excess of the bulk density ρ1. It is further enhanced within the channel confined between two surfaces because a particle of component 1 contacts two surfaces with the result of a more pronounced decrease in the EV for the other particles of component 1. The force is constituted by attractive and repulsive components. They are both very strong and comparable in magnitude, and the net force is much weaker. At h = d1, for example, ρ1,Contact is considerably high within the channel with the result that the repulsive component dominates, thus making the net force repulsive. At h = 0, vanishing of ρ1,Contact occurs for a significant portion of each solute surface, leading to an attractive net force. (See our earlier publications [33,34] for more details.) In the AO theory [31,32], nonzero ρ1,Contact is always set equal to ρ1: the AO force is a purely attractive one with a range of d1. For 0 b h b dS, the remarkable underestimation occurs for both of the attractive and repulsive components. The AO theory benefits from the resultant

Fig. 3. Reduced number density profiles of solvent and cosolvent at contact, GS and GC, respectively, plotted against bulk cosolvent concentration XC = ρC / (ρS + ρC). (a) GS in systems 1 and 2. (b) GC in systems 1 and 2. See Table 1 for the notations.

408

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

  3 3 Γ C ¼ ρC dS GC = ρC dS GC

Pure

:

ð8bÞ

(ρSd3SGS)Pure is ρSd3SGS in the case where the cosolvent is present at infinite dilution (equal to ρSd3S GS in pure solvent): (ρSd3S GS)Pure = 4.184. ρCd3SGC is a measure of the number of cosolvent particles in contact with the solute NContact. (ρCd3S GC)Pure is ρCd3S GC in the case where the solvent is present at infinite dilution (equal to ρCd3S GC in pure cosolvent): (ρCd3S GC)Pure = 0.5008 for dC/dS = 2 and (ρCd3S GC)Pure = 0.0578 for dC/dS = 4. Fig. 4 shows ΓS and ΓC plotted against XC. Intuitively, 0 ≤ ΓS ≤ 1 and 0 ≤ ΓC ≤ 1. Nevertheless, ΓC reaches ~2 at XC ~ 0.01 in system 2. For dC / dS = 4, with a gradual increase in XC, NContact increases, becomes highest at XC ~ 0.01, and turns to a decreasing function beyond this concentration. Higher XC may lead to larger NContact, but higher XS drives the cosolvent particles to contact the solute more strongly, and the latter effect is essential. If the bulk solvent-cosolvent mixture undergoes phase separation, the IET gives no solutions in the unstable region where the mixture cannot exist as a single phase even in the metastable state. In the numerical calculation using the IET, as XC increases from 0 or decreases from 1, the concentration-concentration structure factor at zero wave vector S(0) expressed by Sð0Þ ¼ X S X C ½1 þ ðρS þ ρC ÞX S X C fH SS ð0Þ þ HCC ð0Þ−2H SC ð0Þg

ð9Þ

diverges at the two spinodal values of XC [35]. Here, hij(r) is the total correlation function of the i-j pair: solvent-solvent (SS), cosolventcosolvent (CC) or solvent-cosolvent (SC) pair and Hij(k) is the Fourier transform of hij(r). Hij(0) is equivalent to the Kirkwood-Buff integral [36]. In this study, no such divergent behavior is encountered and the

IET gives converged solutions for all XC even in the case of dC / dS = 4. We find that S(0) is always very close to zero: for instance, S(0) = 0.0088 at XC = 0.012 for dC/dS = 4. The bulk mixture considered in this study is thermodynamically stable, and the anomalous behavior mentioned above originates from the introduction of solute-solvent and solute-cosolvent correlations. 3.5. Examination of thermodynamic consistency of integral equation theory For the bulk solvent-cosolvent mixture without the solute in system 2, we examine the thermodynamic consistency of the IET using Eqs. (1a) and (1b) as the bridge function. The dimensionless partial inverse and χ*−1 compressibilities, χ*−1 S C , respectively, can be expressed as [37] n    o 3 3 3 3 ρS dS C SS ð0Þ=dS þ ρC dS C SC ð0Þ=dS ;

ð10aÞ

n    o 3 3 3 3 χC −1 ¼ 1− ρS dS C SC ð0Þ=dS þ ρC dS C CC ð0Þ=dS :

ð10bÞ

χS −1 ¼ 1−

Here, cij(r) is the total correlation function of the i-j pair: solvent-solvent (SS), cosolvent-cosolvent (CC) or solvent-cosolvent (SC) pair and Cij(k) is the Fourier transform of cij(r). On the other hand, the dimensionless partial inverse compressibilities can also be obtained through numerical differentiation of the virial pressure P with respect to ρS and −1 ρC (they are denoted by (χ*−1 S )Virial and (χ*C )Virial), respectively [37]: P  ¼ dS P=ðkB T Þ 3

 2 3 3 3 g SS ðdS Þ ¼ ρS dS þ ρC dS þ ð2π=3Þ ρS dS    3 3 3 þ2 ρS dS ρC dS ðdSC =dS Þ g SC ðdSC Þ  2 3 3 ðdC =dS Þ g CC dC ; þ ρC dS

f

ð11Þ

g

 

χS −1 χC −1



¼ ∂P  =∂ρS  ; ρS  ¼ ρS dS ;

ð12aÞ

¼ ∂P  =∂ρC  ; ρC  ¼ ρC dS :

ð12bÞ

3

Virial

 Virial

3

Here, dSC = (dS + dC) / 2, the partial derivatives in Eqs. (12a) and (12b) are performed with constant ρC* and T and with constant ρS* and χ*−1 must be and T, respectively, and gij(r) = hij(r) + 1. χ*−1 S C −1 −1 equal to (χ*S )Virial and (χ*C )Virial, respectively. However, since the closure equation is not exact, the equality does not hold, which is referred to as “thermodynamic inconsistency”. In Table 4, under some representative conditions, we present the −1 −1 −1 values of χ*−1 S , χ*C , (χ*S )Virial, and (χ*C )Virial. The condition of XC = 0.012 in the case of dC/dS = 4 is included. As observed in the table, χ*−1 S −1 −1 and χ*−1 C are almost equal to (χ*S )Virial and (χ*C )Virial, respectively. In usual cases, the closure equation includes a parameter or parameters which are adjusted so that the thermodynamic consistency can be satisfied in terms of the total compressibility or partial compressibilities [20, 21,37]. By contrast, the closure equation incorporating Eqs. (1a) and Table 4 Dimensionless partial inverse compressibilities obtained via two different routes using the bridge function expressed by Eqs. (1a) and (1b).

Fig. 4. ΓS and ΓC (see Eqs. (8a) and (8b)) plotted against bulk cosolvent concentration XC = ρC / (ρS + ρC). (a) ΓS in systems 1 and 2. ΓS → 1 as XC → 0 (1 − XC = XS → 1). (b) ΓC in systems 1 and 2. ΓC → 1 as XC → 1. See Table 1 for the notations.

ρSd3S

ρCd3S

χ*−1 S

χ*−1 C

(χ*−1 S )Virial

(χ*−1 C )Virial

0.6677000 0.5397000 0.4117000 0.2837000 0.1557000 0.0091400

0.0010000 0.0030000 0.0050000 0.0070000 0.0090000 0.0112900

17.9 14.0 10.6 7.67 5.27 3.08

542 388 262 162 85.4 23.2

17.8 13.9 10.5 7.68 5.29 3.10

543 390 264 163 86.0 23.1

The bulk solvent-cosolvent mixture without the solute in system 2 is considered (dC/dS = and χ*−1 4). χ*−1 S C are the dimensionless partial inverse compressibilities calculated via −1 Eqs. (10a) and (10b), and (χ*−1 S )Virial and (χ*C )Virial are those calculated via Eqs. (11), (12a) and (12b).

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

(1b) as the bridge function includes no such adjusted parameters. Nevertheless, the IET using Eqs. (1a) and (1b) is almost perfectly consistent. For −1 −1 −1 comparison, the values of χ*−1 S , χ*C , (χ*S )Virial, and (χ*C )Virial obtained from the HNC are presented in Table 5: it is indicated that the HNC suffers significant inconsistency. 3.6. Relevance to real systems The structure and properties of a dense fluid are determined primarily by the repulsive part of the interaction potential [17]. This is why a hard-sphere system has been investigated as the most fundamental one in liquid state theory [17]. Of course, in the real system, vdW and electrostatic interactions are always present. However, the aim of our study is to examine purely entropic effects by shutting off these interactions and show that the effects are quite large for a mixture of particles with high size asymmetry. We believe that the qualitative aspects of the results are not vitiated by the incorporation of these attractive interactions. The reason why the counterintuitive behavior revealed in this study has not been reported in theoretical and experimental studies should be the following: the cosolvent concentration at a surface has not been analyzed or measured by systematically changing the cosolvent concentration in the bulk for a mixture of particles with sufficiently high size asymmetry. A system which can be modeled as hard spheres does exist in colloidal suspensions. Let us consider colloidal particles with highly charged surfaces (e.g., polystyrene microspheres) immersed in water. When NaCl is added to water in a concentration of ~0.01 mol/l, the electrostatic repulsion between particles is screened and the resultant potential can be approximated by a hard-sphere one: the effective diameter of a particle is the sum of the diameter of the particle itself and the Debye length. The aqueous solution can be viewed as inert background. For a similar reason, the interaction between a glass wall and a colloidal particle can be regarded as the hard-body one. Such colloidal systems have been investigated rather extensively in experiments [38–40]. A colloidal suspension containing particles of significantly different sizes near a glass wall is a good model for testing the counterintuitive result revealed in this study. 4. Conclusion In this study, we consider a very large solute immersed in a solventcosolvent binary mixture at infinite dilution. The solute, cosolvent, and solvent molecules are modeled as neutral hard spheres with diameters dU, dC, and dS (dU / dS = 20; dC / dS = 2 or 4), respectively, interacting through simple hard-sphere potentials with no attractive tails, so that the system behavior is purely entropic in its origin. The bulk binary mixture is thermodynamically stable as a single phase in the whole cosolvent-concentration range: it undergoes no phase separation. The number density profiles of solvent and cosolvent near the solute are calculated by an integral equation theory [17] in which a reliable bridge function [25] is incorporated. The bridge function was especially developed for systems with high size asymmetry. Many different values of XC, Table 5 Dimensionless partial inverse compressibilities obtained via two different routes using the HNC. ρSd3S

ρCd3S

χ*−1 S

χ*−1 C

(χ*−1 S )Virial

(χ*−1 C )Virial

0.6677000 0.5397000 0.4117000 0.2837000 0.1557000 0.0091400

0.0010000 0.0030000 0.0050000 0.0070000 0.0090000 0.0112900

13.6 10.8 8.21 5.97 4.10 2.53

410 297 203 126 66.5 17.6

23.3 17.9 13.3 9.45 6.29 3.49

748 532 355 216 112 30.3

The bulk solvent-cosolvent mixture without the solute in system 2 is considered (dC/dS = and χ*−1 4). χ*−1 S C are the dimensionless partial inverse compressibilities calculated via −1 Eqs. (10a) and (10b), and (χ*−1 S )Virial and (χ*C )Virial are those calculated via Eqs. (11), (12a) and (12b).

409

the cosolvent mole fraction in the bulk, are tested. We note that in this study the component with a larger particle diameter is categorized as the cosolvent and cases of XS b XC (XS = 1 − XC) are also analyzed to investigate the basic behavior of the entropic excluded-volume effect. The solvent and cosolvent number densities at the solute surface become higher than in the bulk: the solvent and the cosolvent are enriched near the solute. Some of the important results are recapitulated below. (1) When the cosolvent is added to the solvent at very low concentration, the cosolvent is significantly enriched near the solute: the degree of enrichment is much higher than in pure cosolvent, and this trend becomes stronger for larger dC. When the solvent is added to the cosolvent at very low concentration, the solvent is not significantly enriched near the solute: the degree of enrichment is lower than in pure solvent, and this trend becomes stronger for larger dC. (2) The Asakura-Oosawa (AO) theory [31,32], which is routinely applied to the estimation of the solute-solute entropic potential at contact, gives a fairly accurate result in a quantitative sense only when the diameters of two solute particles are not significantly different. The AO estimation becomes less accurate as the size asymmetry of two solute particles increases. (3) As XC gradually increases, the number of cosolvent particles in contact with the solute NContact becomes larger. When dC/dS is sufficiently large, however, NContact becomes largest at a surprisingly low value of XC: beyond this value of XC, NContact continues to decrease. For dU/dS = 20, dC/dS = 4, XC ~ 0.01 and the maximum value of NContact is approximately twice larger than NContact in pure cosolvent. This anomalous behavior does not arise from the solvent-cosolvent correlations because the bulk mixture is always thermodynamically stable. It is ascribed to the introduction of solute-solvent and solute-cosolvent correlations. (4) It has been corroborated that the integral equation theory using the bridge function expressed by Eqs. (1a) and (1b) satisfies the thermodynamic consistency almost perfectly for the bulk solvent-cosolvent mixture with dC/dS = 4. Of course, the results obtained in this study are relevant not only to protein-water-cosolvent systems but also to general fluid mixtures at surfaces. In the real system, the degree of cosolvent enrichment near the solute is influenced not only by the entropic effect analyzed in this study but also by energetic factors such as cosolvent-solvent and cosolvent-solute interaction potentials. What we wish to emphasize is the following: the entropic effect can be remarkably large; and without an analysis based on statistical mechanics it is difficult to quantitatively conjecture the dependence of the degree of cosolvent enrichment on dU/dS, dC/dS, and XC. There can be three subjects to be pursued further. First, a more extensive parametric study considering additional combinations of (dU/dS, dC/dS) (e.g., (dU/dS, dC/dS) = (20, 3) and (10, 4)) is strongly desired. Second, the contributions from van der Waals (vdW) and electrostatic interaction potentials can be investigated by accounting for only vdW interaction potential in the first step and both of the interaction potentials in the second one in the model on a step-by-step basis. Third, it is interesting to calculate the potential of mean force (PMF) between very large solutes, especially under the condition: dU/dS = 20, dC/dS = 4, XC ~ 0.01. Acknowledgment This work was supported by Grant-in-Aid for Scientific Research (B) (No. 17H03663) from Japan Society for the Promotion of Science (JPSJ) to M. K. References [1] P.H. von Hippel, K.-Y. Wong, On the conformational stability of globular proteins, J. Biol. Chem. 240 (1965) 3909–3923.

410

M. Kinoshita, T. Hayashi / Journal of Molecular Liquids 247 (2017) 403–410

[2] T.T. Herskovits, H. Jaillet, B. Gadegbeku, On the structural stability and solvent denaturation of proteins. II. Denaturation by the ureas, J. Biol. Chem. 245 (1970) 4544–4550. [3] N. Poklar, G. Vesnaver, S. Lapanje, Thermodynamics of denaturation of α-chymotrypsinogen A with alkylureas, J. Protein Chem. 14 (1995) 709–719. [4] G. Graziano, How does trimethylamine N-oxide counteract the denaturing activity of urea? Phys. Chem. Chem. Phys. 13 (2011) 17689–17695. [5] R. Sarma, S. Paul, Exploring the molecular mechanism of trimethylamine-N-oxide's ability to counteract the protein denaturing effects of urea, J. Phys. Chem. B 117 (2013) 5691–5704. [6] Y. Yamamori, R. Ishizuka, Y. Karino, S. Sakuraba, N. Matubayasi, Interaction-component analysis of the hydration and urea effects on cytochrome c, J. Chem. Phys. 144 (2016) 085102. [7] S. Murakami, M. Kinoshita, Effects of monohydric alcohols and polyols on the thermal stability of a protein, J. Chem. Phys. 144 (2016) 125105. [8] S. Murakami, T. Hayashi, M. Kinoshita, Effects of salt or cosolvent addition on solubility of a hydrophobic solute in water: relevance to those on thermal stability of a protein, J. Chem. Phys. 146 (2017) 055102. [9] Y. Harano, M. Kinoshita, Translational-entropy gain of solvent upon protein folding, Biophys. J. 89 (2005) 2701–2710. [10] M. Kinoshita, A new theoretical approach to biological self-assembly, Biophys. Rev. 5 (2013) 283–293. [11] E. López-Sánchez, C.D. Estrada-Álvarez, G. Pérez-Ángel, J.M. Méndez-Alcaraz, P. González-Mozuelos, R. Castañeda-Priego, Demixing transition, structure, and depletion forces in binary mixtures of hardspheres: the role of bridge functions, J. Chem. Phys. 139 (2013) 104908. [12] A. Suematsu, A. Yoshimori, R. Akiyama, Effects of interactions between depletants in phase diagrams of binary hard-sphere systems, EPL 116 (2016) 38004. [13] T. Biben, P. Bladon, D. Frenkel, Depletion effects in binary hard-sphere fluids, J. Phys. Condens. Matter 8 (1996) 10799–10821. [14] R. Dickman, P. Attard, V. Simonian, Entropic forces in binary hard sphere mixtures: theory and simulation, J. Chem. Phys. 107 (1997) 205–213. [15] M. Kinoshita, S. Iba, K. Kuwamoto, M. Harada, Interaction between macroparticles in Lennard-Jones fluids or in hard-sphere mixtures, J. Chem. Phys. 105 (1996) 7177–7183. [16] R. Roth, M. Kinoshita, Depletion potential between large spheres immersed in a multi-component mixture of small spheres, J. Chem. Phys. 125 (2006) 084910. [17] J.-P. Hansen, L.R. McDonald, Theory of simple liquids, 3rd ed. Academic, London, 2006. [18] P. Attard, G.N. Patey, Hypernetted-chain closure with bridge diagrams. Asymmetric hard sphere mixtures, J. Chem. Phys. 92 (1990) 4970–4982. [19] L. Verlet, Integral equations for classical fluids. I. The hard-sphere case, Mol. Phys. 41 (1980) 183–190. [20] F.J. Rogers, D.A. Young, New, thermodynamically consistent, integral equation for simple fluids, Phys. Rev. A 30 (1984) 999–1007. [21] G. Zerah, J.-P. Hansen, Self-consistent integral equations for fluid pair distribution functions: another attempt, J. Chem. Phys. 84 (1986) 2336–2343.

[22] A. Malijevský, S. Labík, The bridge function for hard spheres, Mol. Phys. 60 (1987) 663–669. [23] D.-M. Duh, D. Henderson, Integral equation theory for Lennard-Jones fluids: the bridge function and applications to pure fluids and mixtures, J. Chem. Phys. 104 (1996) 6742–6754. [24] N. Choudhury, S.K. Ghosh, Integral equation theory of penetrable sphere fluids: a modified Verlet bridge function approach, J. Chem. Phys. 119 (2003) 4827–4832. [25] M. Kinoshita, Interaction between surfaces with solvophobicity or solvophilicity immersed in solvent: effects due to addition of solvophobic or solvophilic solute, J. Chem. Phys. 118 (2003) 8969–8981. [26] A. Santos, R. Fantoni, A. Giacometti, Thermodynamic consistency of energy and virial routes: an exact proof within the linearized Debye–Hückel theory, J. Chem. Phys. 131 (2009) 181105. [27] T. Morita, Theory of classical fluids: hyper-netted chain approximation. IIIa. A new integral equation for the pair distribution function, Prog. Theor. Phys. 23 (1960) 829–845. [28] T. Morita, K. Hiroike, A new approach to the theory of classical fluids. III: general treatment of classical systems, Prog. Theor. Phys. 25 (1961) 537–578. [29] M. Kinoshita, M. Harada, Numerical solution of the RHNC theory for fluids of nonspherical particles near a uniform planar wall. A study with dipolar hard spheres as an example system, Mol. Phys. 79 (1993) 145–167. [30] M. Kinoshita, D.R. Bérard, Analysis of the bulk and surface-induced structure of electrolyte solutions using integral equation theories, J. Comput. Phys. 124 (1996) 230–241. [31] S. Asakura, F. Oosawa, On interaction between two bodies immersed in a solution of macromolecules, J. Chem. Phys. 22 (1954) 1255. [32] S. Asakura, F. Oosawa, Interaction between particles suspended in solutions of macromolecules, J. Polym. Sci. 33 (1958) 183–192. [33] M. Kinoshita, Spatial distribution of a depletion potential between a big solute of arbitrary geometry and a big sphere immersed in small spheres, J. Chem. Phys. 116 (2002) 3493–3501. [34] M. Kinoshita, Roles of entropic excluded-volume effects in colloidal and biological systems: analyses using the three-dimensional integral equation theory, Chem. Eng. Sci. 61 (2006) 2150–2160. [35] A.B. Bhatia, D.E. Thornton, Structural aspects of the electrical resistivity of binary alloys, Phys. Rev. B 2 (1970) 3004–3012. [36] J.G. Kirkwood, F.P. Buff, The statistical mechanical theory of solutions. I, J. Chem. Phys. 19 (1951) 774–777. [37] T. Biben, J.-P. Hansen, Phase separation of asymmetric binary hard-sphere fluids, Phys. Rev. Lett. 66 (1991) 2215–2218. [38] P.D. Kaplan, J.L. Rouke, A.G. Yodh, D.J. Pine, Entropically driven surface phase separation in binary colloidal mixtures, Phys. Rev. Lett. 72 (1994) 582–585. [39] A.D. Dinsmore, A.G. Yodh, D.J. Pine, Phase diagrams of nearly hard-sphere binary colloids, Phys. Rev. E 52 (1995) 4045–4057. [40] A.D. Dinsmore, A.G. Yodh, D.J. Pine, Entropic control of particle motion using passive surface microstructures, Nature 383 (1996) 239–242.