Analysis of Molecular Aggregation Equilibria Using ... - ACS Publications

Article pubs.acs.org/JPCB

Analysis of Molecular Aggregation Equilibria Using Random Mixing Statistics Blake M. Rankin and Dor Ben-Amotz* Department of Chemistry, Purdue University, 560 Oval Drive, West Lafayette, Indiana 47906, United States ABSTRACT: Aggregation processes in both the gas phase and aqueous solutions are analyzed by comparing aggregate size distributions obtained from molecular dynamics simulations with analytical predictions pertaining to a nonaggregating random mixture. The latter predictions are obtained by using the binomial distribution to predict the statistical properties of a uniformly mixed solution containing molecules of the same size and concentration as those in the solution of interest. Simulations are performed on systems containing neopentane dissolved in either methane, aqueous methanol, or aqueous NaI solutions. Comparisons of the theoretical and simulation results are used to both classify and quantify the influence of intermolecular interactions on such aggregation processes, including the equilibrium constants and thermodynamic functions pertaining to the partitioning of molecules between the bulk and first coordination shell of neopentane. Although the present results are primarily intended to describe and illustrate the random mixing analysis strategy, they also imply that neopentane has a greater tendency to aggregate with methane in the dense gas phase than with either methanol or iodide ions in aqueous solutions.

■

INTRODUCTION Molecular interactions of biological and environmental interest are often driven by a delicate balance of noncovalent forces and thermal fluctuations. These include, for example, processes involving hydrophobic aggregation and specific ion (Hofmeister) interactions. Such processes typically occur at relatively high concentrations, at which one would expect there to be a substantial number of intermolecular contacts even in a randomly mixed (nonaggregating) system. Thus, in order to determine if a particular process involves the enhanced formation or breakup of aggregates, it is necessary to understand the aggregation statistics associated with the corresponding randomly mixed reference system. We have recently suggested a means of doing that, by making use of the binomial distribution to predict cluster size distributions in random mixtures resembling aqueous t-butyl alcohol solutions.1 Here, we generalize this strategy to other types of aggregation processes and compare the resulting predictions with molecular dynamics (MD) simulation results pertaining to both dense gas mixtures and aqueous solutions containing neopentane, methane, methanol, and/or NaI. Our results serve to validate the random mixing model and show how this strategy can be used to both classify and quantify molecular aggregation equilibria. This work is motivated in part by recent experiments and simulations suggesting that the degree of aggregation in some aqueous alcohol solutions may not exceed that in a random mixture, and under some conditions, there may be even fewer direct hydrophobic contacts than in a random mixture.1,2 The work is also motivated by recent debates regarding whether ions of various types accumulate or are expelled from molecular interfaces, such as the first hydration shells of hydrophobic groups dissolved in an aqueous salt solution.3−9 In all such systems, a careful analysis is required in order to distinguish © XXXX American Chemical Society

small deviations from random mixing and thus determine whether a particular alcohol or ion does or does not have an affinity to coordinate with a particular solute or subgroup. We define a randomly mixed solution as one in which the concentration of each component is uniformly distributed through the system. In other words, all of the pair distribution functions in such a system are assumed to be equal to 1, outside of the volume excluded by the core of each molecule. Under such conditions, the bulk concentration of each molecule is all that is required in order to predict the probability of finding that molecule within a given region. For example, in a random mixture with a solute concentration of [c], there will be an average of exactly V[c] solute molecules with any volume V. If the volume of interest is smaller than the solute’s partial molar volume V ≤ V̅ , then V[c] becomes the probability of finding a single solute molecule within that small volume. Note that V̅ [c] = 1 in a single-component fluid because in this case, V̅ = 1/[c], and therefore, in any solution of uniform concentration, V̅ [c] < 1 for each component. Our current interest is in establishing the probability that the first coordination shell around a given solute molecule will contain a particular number of other molecules, which we will from here on refer to as “ligands”. The maximum number n of ligands that can be accommodated within the first coordination shell of volume Vshell should be close to that in a close-packed structure, in which the solute is entirely surrounded by ligands. In a random mixture, the average number of ligands in the coordination shell is exactly Vshell[c]. One may extend this line Special Issue: Michael D. Fayer Festschrift Received: June 28, 2013 Revised: August 16, 2013

A

dx.doi.org/10.1021/jp406413b | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

of reasoning by making use of the binomial distribution to predict the probabilities that there will be a particular number of ligands within the first coordination sphere around a given solute in a random mixture (as further described in the next section). By doing so, one may establish the random mixing aggregate size distribution as a benchmark with respect to which the aggregate size distribution in a real mixture can be classified and quantified. In other words, we may use this strategy not only to determine how the aggregate size distributions observed in a given experiment or computer simulation differ from random mixing predictions but also to quantify the increase or decrease in the local ligand concentration, relative to that in the surrounding solution. Moreover, the latter concentration may be used to determine the partition coefficient of the ligand between the bulk solution and the solute coordination shell and thus to obtain the corresponding thermodynamic functions.

The above expression is quite general as it may be used to predict random mixing probability distributions within any volume of interest for solutions containing any number of components. For all of the solute−ligand aggregation processes described in this work, we equate n with the maximum number of ligands that can close-pack around the solute. The following geometric argument may be used to obtain a physically reasonable estimate of n from the ratio of the solute−ligand accessible surface area ASA12 to the surface area occupied by a single ligand A2. ⎞2 ⎛ σ1 ASA12 n≈ ≈ π ⎜ + 1⎟ A2 ⎠ ⎝ σ2

The second approximate equality is obtained when the solute and ligand are described as spheres with effective diameters of σ1 and σ2, respectively, in which case ASA12 ≈ π(σ1 + σ2)2 and A2 ≈ σ22 approximate the portion of the accessible surface area that is occupied by a single bound ligand. When eq 2 produces noninteger values, the maximum coordination number n can be approximated by rounding down to the nearest integer. Note that eq 2 is consistent with the exact close-packed values of n = 12 when σ1 = σ2 and n = 3 when σ1 = 0 (because no more than three spheres can close-pack around a point in threedimensional space). Moreover, eq 2 predicts that n > 12 when σ1 > σ2, which is consistent with the fact that more than 12 small spheres can close-pack around a larger sphere. It is also noteworthy that replacing π by 3 in eq 2 produces exact agreement with the above close-packing values. As previously discussed, the random mixing probability p of finding a ligand of concentration [c2] within a subcell of volume (1/n)Vshell is necessarily given by the first of the following equalities

■

RANDOM MIXING MODEL The random mixing approximation may in general be utilized to predict the probabilities of finding any number of components (of either the same or different molecular structure) in any volume V within the system, given only the concentration of each component. The volume of interest may also, in general, be partitioned into n subcells of equal volume V/n. The probability of finding a component whose bulk concentration is [ci] in each subcell is p = (1/n)V[ci]. Because p is inversely proportional to n, the precise value of n has little influence on the predicted aggregate size distribution, except when n is so small that it is comparable to the maximum number of molecules that can fit within the volume V (which, as previously noted, is the smallest physically reasonable value of n). Any value of n larger than or equal to this value will ensure both that there is a negligible probability that more than one molecule will be found within each subcell and that the predicted aggregate size distribution is relatively insensitive to n. When applied to the sorts of solute−ligand aggregation processes described in this work, the above generalized formulation may be implemented by taking n and VShell to be the fundamental parameters of the random mixing aggregation process, and thus, p = (1/n)VShell[c2], where [c2] is the bulk concentration of the ligand molecules. Choosing n to be equal to the maximum ligand coordination number (as further described below) leads to random mixing aggregate size predictions that reflects the steric constraint on the coordination number, while larger values of n produce quite similar random mixing predictions pertaining to a situation in which there is no absolute upper bound placed on the number of ligands that can occupy the solute’s coordination shell. However, because the concentration of any molecule also cannot exceed an upper bound imposed by its finite core volume, the n-independent predictions obtained in the large n limit are themselves necessarily consistent with steric closepacking constraints. The probability P(k) that k subcells will be occupied by a ligand in a random mixture may be determined using the following binomial distribution expression. P(k) =

n! pk (1 − p)n − k k ! (n − k )!

(2)

1 1 VShell[c 2] = [c 2] n n 4π 3 3 (r2 − r1 )[c 2] ≈ 3n

p=

∫0

r2

g0(r )4πr 2 dr (3)

Note that when eq 2 is used to obtain n, then (1/n)Vshell is typically similar to (but slightly smaller than) the partial molar volume of each ligand V̅ 2, and thus, p = (1/n)VShell[c2] ≤ V̅ 2[c2] ≤ 1. The second equality in eq 3 is obtained by equating the solute−ligand pair distribution function g0(r) (pertaining to the central atoms of the solute and ligand molecules) with the random mixing pair distribution function that is equal to 1 outside of the solute−ligand contact core volume. Because realistic interaction potentials have soft-core volumes, one may use the following functional form for g0(r) to approximate the soft leading edge of such solute−ligand pair distribution functions. ⎡ ⎛ r ⎞m ⎤ g0(r ) = exp⎢ −ln(2)⎜ 1 ⎟ ⎥ ⎝r⎠ ⎦ ⎣

(4)

The factor of ln(2) assures that g0(r) = 1/2 when r = r1. In this work, we have assigned the inverse power law exponent to m = 18 as this exponent approximates the steepness of the repulsive core of the Lennard-Jones potential of argon10 and also reasonably represents the leading edge of the g(r) simulation results obtained in this work. Increasing the value of m produces a steeper g0(r) that approaches that of hard spheres with a contact distance of r1. The last expression in eq 3

(1) B



■

RESULTS AND DISCUSSION Neopentane in Methane. Figure 1 shows both MD simulation results and random mixing predictions pertaining to a nonideal gas system containing neopentane dissolved in methane, at both 293 (upper panels) and 1000 K (lower panels). The methane concentration range of 0.5 ≤ [c2] ≤ 15 M spans that of a nearly ideal gas to dense fluid as 15 M is 50% above the critical density of methane (and about half of its triple point density). Note that the aggregate size distribution functions on the right-hand side of Figure 1 contain more information than the pair distribution functions on the lefthand side of Figure 1. This is because g(r) only provides information about the average methane concentration at various distances from neopentane, while P(k) reveals not only the average but also the width (and shape) of the corresponding C(CH3)4·[CH4]k aggregate size distributions, within a sphere of radius r2 around the central carbon of neopentane, which specifies the extent of the first coordination shell. The value of r2 may reasonably be assigned a value of 7 Å, as that is the location of the first minimum in the neopentane− methane g(r) at a methane concentration of 15 M. The aggregate size distributions in panels B and D of Figure 1 show the probabilities that exactly k methane molecules will occupy the first coordination shell of neopentane. The distributions indicated by open points (and dashed lines) represent random mixing predictions obtained using eqs 1−3. The latter distributions are obtained assuming either a soft- or hard-core random mixing g0(r), as indicated by the dashed and dotted curves in panels A and C. The aggregate size distributions indicated by solid points (and solid curves) in panels B and D are the corresponding MD simulation results (again obtained using the same coordination shell radius of r2 = 7 Å, as further discussed below). The pair distribution functions g(r) shown in panels A and B of Figure 1 clearly depend on both the concentration of methane and the temperature of the system. At low methane concentration, g(r) necessarily approaches e−V(r)/kBT, where V(r) is the spherically averaged neopentane−methane pair potential (where kB is Boltzmann’s constant and T is the absolute temperature).22 At T = 293 K and [c2] = 15 M, g(r) shows evidence of the formation of both a first and second coordination shell of methane around neopentane. The value of r1 = 4.4 Å obtained from the leading edge of g(r) is close to that obtained using r1 ≈ (1/2)(σ1 + σ2) ≈ 4.54 Å (where σ1 = 5.50 Å and σ2 = 3.58 Å are the neopentane and methane diameters, respectively, obtained from the compressibility of the corresponding pure fluids).11 The value r2 = 7.0 Å obtained from the first minimum in the neopentane−methane g(r) at 15 M and 293 K is somewhat larger than r2 ≈ (σ1/2) + σ2 ≈ 6.33 Å, which may be viewed as a lower bound to the extent of the first coordination shell of methane around neopentane. Other values of r1 and r2, ranging from anywhere within the leading edge of g(r) to somewhere near the first minimum in g(r), may also reasonably be used to define the inner and outer boundaries of the first coordination shell (as further discussed below in the subsection entitled Robustness of the Random Mixing Analysis). The aggregate size distributions derived from the MD simulations, shown in panels B and D of Figure 1, are obtained by counting the number of times that exactly k methane molecules are found within a sphere of radius of r2 = 7 Å about the central carbon of neopentane. The histograms in Figure 1

pertains to such a hard-sphere g0(r) and leads to random mixing predictions that are virtually indistinguishable from those obtained using the soft-core g0(r) integral in eq 3, as long as r1 is taken to realistically represent the solute−ligand contact distance. The latter contact distances can be obtained from the leading edge of the solute−ligand g(r) or, when that is not feasible, from the mean value of the effective hard-sphere (or van der Waals) diameters of the solute σ1 and ligand σ2 molecules, r1 ≈ (1/2)(σ1 + σ2).11−13 Our results indicate that physically reasonable variations of both r1 and r2 have a relatively small influence on the quantitative random mixing predictions and MD simulations (as further discussed below in the subsection entitled Robustness of the Random Mixing Analysis). The probabilities given by eq 1 may further be used to obtain the corresponding excess chemical potential of a cluster that contains exactly k ligands. μk = −kBT ln[P(k)]

Article

(5)

This excess chemical potential is that of an aggregate in a randomly mixed solution, measured relative to the absolute chemical potentials of the solute (of concentration [c1]) and k ligands (of concentration [c2]) in the same solution. This excess chemical potential is typically a nonlinear function of k, with a minimum at the most probable aggregate size, which is in turn close to the average aggregate size ⟨k⟩ = pn = (1/n)VShell[c2]n = VShell[c2]. The curvature of μk about the minimum reflects the entropic cost associated with fluctuations that vary with the local ligand concentration in the coordination shell away from the surrounding bulk ligand concentration, in a randomly mixed solution.

■

SIMULATION METHODS MD simulations were performed using GROMACS on dense gas and aqueous systems.14 The gas-phase system contained 50 methane molecules and one neopentane C(CH3)4 molecule. Methane CH4 concentrations ranging from 0.5 to 15 M were obtained by sequentially decreasing the box size. The aqueous methanol simulations consisted of a single neopentane molecule and methanol dissolved in approximately 500 TIP4P water molecules, while the aqueous NaI simulations contained approximately 1000 TIP4P water molecules15 and a sufficient number of either methanol CH3OH molecules or sodium Na+ and iodide I− ions to produce a concentration of ∼5 M. Both neopentane and methane were represented using OPLS-AA potentials,16 and all bond lengths and angles were constrained using the SHAKE algorithm17 to geometries optimized using B3LYP-6-31G.18 A velocity-rescale Berendsen thermostat19 and Berendsen barostat20 were used to perform NPT simulations at 293 K and 1 bar for the aqueous systems and at either 293 or 1000 K for the gas-phase system (pressures were kept constant to achieve the desired concentrations). Thermodynamic results in Figure 4 were obtained from simulations performed at T = 273, 293, and 313 K. The particle mesh Ewald (PME) method21 was used for electrostatic interactions, and a cutoff of less than half of the box length was used for Lennard-Jones interactions. Pair distribution functions were calculated from the central carbon atom of neopentane to the carbon atom of methane or methanol and to the iodide ion. All systems were pre-equilibrated for at least 200 ps. Results were obtained from MD simulations of 10 ns duration. C



Article

Figure 1. Comparison of random mixing predictions and MD simulations for a system consisting of neopentane dissolved in methane at various methane concentrations of 0.5 ≤ [c2] ≤ 15 M and temperatures of 293 and 1000 K. The dashed curves (and open points) represent random mixing predictions, and the solid curves (and closed points) are MD simulation results. The dashed and dotted pair distribution functions in panels (A) and (C) represent the corresponding hard- and soft-core g0(r) functions (both of which produce essentially identical random mixing aggregate size distribution predictions).

the left, and the aggregate size distributions are shown on the right. The values of r1 and r2 used to obtain these results (and given in the caption of Figure 2) are again obtained from the leading edge and first minimum of the corresponding g(r). Physically reasonable variations of r1 and r2 do not alter the following conclusions (as further discussed in the subsection entitled Robustness of the Random Mixing Analysis). The results for neopentane in aqueous methanol (upper two panels in Figure 2) differ in several ways from the corresponding results for neopentane dissolved in methane gas (shown in Figure 1). First of all, the pair distribution function in Figure 2A contains both first and second coordination shell peaks, while Figure 1A shows no evidence of a second coordination layer (at the same 5 M ligand concentration). More interestingly, while substantial aggregation of methane around neopentane was found at a methane concentration of 5 M (at 293 K), the aggregate size distributions shown in Figure 2B indicate that there is somewhat less aggregation of methanol around neopentane in water (at the same concentration and temperature). Thus, the MD simulation and random mixing predictions imply the van der Waals interactions between neopentane and methane in the gas phase provide a greater driving force for aggregation than hydrophobic interactions between neopentane and methanol in liquid water. The results obtained for neopentane in 5 M aqueous NaI, which are shown in the bottom two panels of Figure 2, tell a quite different story. In this case, the neopentane−I− g(r) has a significantly smaller maximum value, and the corresponding aggregate size distributions indicate that there are slightly fewer I− ions in the first hydration shell of neopentane than in a random mixture. In other words, these simulations indicate that I− is mildly expelled from the first hydrophobic hydration shell of neopentane. Such results are undoubtedly sensitive to the nature of the interaction potential between I− and neo-

all shift to increasing aggregate sizes with increasing methane concentration. These distributions indicate, for example, that at 0.5 M, it is most probable that neopentane will have no methane molecules (k = 0) in its first coordination shell, while at 15 M, the MD simulation results predict that neopentane will have an approximately Gaussian aggregate size distribution peaked at k ≈ 12.5 with a full width at half-maximum of Δk ≈ 5. Comparisons of the random mixing (dashed) and MD simulation (solid) histograms in Figure 1B indicates that there are significantly more methane molecules in the first coordination shell of neopentane than in a randomly mixed system. The latter conclusion holds at all three methane concentrations and is consistent with the fact that the MD g(r) functions are greater than 1 throughout most of the first coordination shell. Thus, van der Waals interactions evidently lead to significant aggregation of methane around neopentane at 293 K. However, when the temperature of the simulations is increased to 1000 K, the peak in g(r) decreases, and the corresponding aggregate size distributions not only shift downward but also reveal that there is no longer as large of a difference between the MD simulation and random mixing predictions at methane concentrations of [c2] ≤ 5 M, while at 15 M, the enhanced structuring of methane around neopentane again leads to significant deviations from random mixing predictions. These results and conclusions are remarkably insensitive to the precise values of r1 and r2 used to obtain the random mixing and MD simulation aggregate size distribution histograms (as further explained below in the subsection entitled Robustness of the Random Mixing Analysis). Neopentane in Aqueous Alcohol and Salt Solutions. Figure 2 compares random mixing and MD simulation results pertaining to aqueous solutions containing one neopentane molecule dissolved in either 5 M methanol or 5 M NaI at 293 K. The corresponding pair distribution functions are shown on D



Article

Figure 2. Comparison of MD simulation results with random mixing predictions for neopentane dissolved in aqueous solutions containing 5 M methanol (A and B) or 5 M NaI (C and D). The values of r1 and r2 obtained from the leading edge and first minimum of g(r) are r1 = 4.4 Å and r2 = 7.06 Å for neopentane in aqueous methanol and r1 = 5.05 Å and r2 = 7.26 Å for the aggregation of I− around neopentane in aqueous NaI.

pentane.9,23−25 However, the conclusion reached using the present simulation results is qualitatively consistent with recent experimental measurements that also imply that I− is expelled from the hydrophobic hydration shells of other hydrophobic molecules and subgroups dissolved in water.3,4 Robustness of the Random Mixing Analysis. To further test and validate the conclusions reached using the present random mixing analysis, it is important to determine the sensitivity of the results to variations in the input parameters r1 and r2 (that specify the extent of the first coordination shell). The most significant variations are those obtained when changing the value of r2 around the first minimum in g(r) as any value of r1 within the leading edge of g(r) has little or no effect on the resulting random mixing predictions. Thus, in Figure 3, we show how the aggregate size distributions obtained from the MD simulation results and random mixing predictions depend on physically reasonable variations in r2 for systems containing neopentane surrounded by either methane, aqueous methanol, or aqueous I− (all at the same 5 M concentration and 293 K). Panels A and B in Figure 3 show results pertaining to neopentane dissolved in a dense gas of methane. The value of r2= 6.33 Å in Figure 3A is that obtained from the experimentally derived effective hard-sphere diameters σ1 and σ2 of neopentane and methane,11 while r2 = 7 Å corresponds to the location of the first minimum in the neopentane−methane g(r) at 15 M (and 293 K). The corresponding aggregate size distributions shown in Figure 3B reveal that even such substantial variations in r2 only produce relatively modest changes in P(k). Most importantly, when the value of r2 is changed, both the random mixing and MD simulation aggregate size distributions shift in the same direction. Thus, the qualitative conclusion that methane aggregates around neopentane at 293 K is unaffected by the precise values of r1 and r2. Panels C and D in Figure 3 pertain to neopentane dissolved in a 5 M aqueous methanol solution. In this case, the smaller value of r2 is obtained using r2 = r1 + (σ2/2) (where σ2 is the

Figure 3. The robustness of conclusions inferred by comparing MD simulations (solid curves in B, D, and F) and random mixing (dashed curves in B, D, and F) are tested by showing how the corresponding aggregate size distributions vary when choosing different physically reasonable values r1 and r2 to define the extent of the first coordination shell of neopentane. (A) and (B) pertain to neopentane in methane, (B) and (C) pertain to neopentane in aqueous methanol, and (E) and (F) pertain to neopentane in aqueous NaI (each obtained at [c2] = 5 M and T = 293 K). Note that the pair distribution functions in panels A, C, and E are identical to those in Figures 1A, 2A, and 2C, respectively.

experimental effective hard-sphere diameter of methanol),11 and the larger value of r2 is obtained from the first minimum in g(r). The aggregate size distributions obtained from MD simulations and random mixing predictions are again found to shift in the same direction when r2 is varied. Comparisons of the results in panels B and D of Figure 3 confirm that the E



Article

aggregation of neopentane with methanol (in water) is invariably slightly less pronounced than that with methane (in the gas phase). Panels E and F in Figure 3 pertain to neopentane dissolved in an aqueous NaI solution. In this case, the smaller and larger values of r2 are obtained from the distance at which g(r) crosses 1 after the first maximum and the location of the first minimum in g(r), respectively. The resulting aggregate size distributions indicate that there is so little difference between the aggregates size distributions obtained from the MD simulation and those from random mixing predictions that reasonable variations in the selected values of r1 and r2 may imply that there is either slight expulsion or slight aggregation of I− in the first hydration shell of neopentane in water. However, these results weigh slightly in favor of expulsion because there is only a very small region over which g(r) exceeds 1, and therefore, most physically reasonable choices of the values of r1 and r2 lead to the conclusion that I− is slightly expelled from first hydration shell of neopentane (according to the present purely classical, nonpolarizable simulations at ambient pressure and temperature). Aggregation Thermodynamics. The above results indicate that the aggregate size distributions obtained from MD simulations may deviate to a greater or lesser extent from the corresponding aggregate size distributions in a random mixture. The latter deviations may be used to quantify the difference between the actual local ligand concentration within the solute’s first coordination shell and that in the corresponding nonaggregating randomly mixed reference system. More specifically, the actual average number of ligands in the first coordination shell ⟨k⟩Shell may be used to determine the effective local ligand concentration within the first coordination shell [c2]Shell as follows. [c 2]Shell

⟨k⟩Shell = VShell

Figure 4. The two left-hand panels show comparisons of aggregate size distributions obtained from MD simulations (points) with those predicted in a random mixture with the same aggregate concentration [c2] = 5 M (open points and dashed lines) and the best fit (solid lines) to the random mixture distribution obtained when treating the ligand concentration as a fit parameter to obtain the local ligand concentration in the first coordination shell of the solute [c2]Shell. The r2 values of neopentane in aqueous methanol and iodide were 7.06 and 7.26 A, respectively (obtained from the first minimum in the pair distribution functions shown in Figure 2). The two right-hand panels show the resulting ΔG/T plotted versus 1/T, whose slopes correspond to ΔH. The resulting values of ΔG, ΔH, and TΔS are shown, along with error bars obtained from five sequential 2 ns simulations.

[c2]Shell is approximate, it can be shown that this procedure becomes exact for a lattice model of such aggregation processes when ligand−ligand interactions are negligible in comparison with solute−ligand interactions (as will be described in a forthcoming publication by the current authors and Ben Widom of Cornell University). Moreover, the resulting expressions for K and ΔG are also exact for a process corresponding to transfer of a ligand from a concentration [c2] to a concentration of [c2]Shell, where ΔG is equivalent to the difference between the excess chemical potentials of the ligand at the two concentrations (relative to that in noninteracting ideal gas systems of the same concentration). The enthalpy ΔH and entropy ΔS associated with the above aggregation process may be obtained from the temperature dependence of ΔG using the following standard thermodynamic relations.

(6)

Alternatively, one may fit the actual aggregate size distribution to that in a random mixture by treating [c2]Shell as a fitting parameter, as illustrated in the left-hand panels in Figure 4. Note that both of these procedures for obtaining an empirical value of [c2]Shell are equivalent to assigning p = (1/n) VShell[c2]Shell to the empirical probability that a ligand will occupy an arbitrary site within the coordination shell. If the actual aggregate size distributions were identical to that in a random mixture of concentration [c2], then eq 6 would yield [c2]Shell = [c2] because ⟨k⟩ = VShell[c2] in a random mixture. More generally, the ratio of [c2]Shell to [c2] is a measure of the partition coefficient of the ligand between the bulk and the solute coordination shell.5 [c ] K = 2 Shell = e−ΔG / RT [c 2]

⎡ ∂(ΔG /T ) ⎤ ΔH = ⎢ ⎥ ⎣ ∂(1/T ) ⎦ P

(7)

and

⎛ ∂ΔG ⎞ ⎟ ΔS = −⎜ ⎝ ∂T ⎠ P

(8)

The right-hand panels in Figure 4 contain ΔG, ΔH, and TΔS results at 293 K obtained from simulations of neopentane dissolved in either aqueous methanol or aqueous I− (each at the same 5 M concentration). The derivatives in eq 8 were determined from linear fits to simulation results obtained at 273, 293, and 313 K. The values of ΔG, ΔH, and TΔS are those obtained from 10 ns simulations. The error bars represent the standard deviation of the means obtained from five sequential 2 ns simulations (extracted from a single 10 ns simulation). The thermodynamic results pertaining to the aggregation of neopentane and methanol in water, shown in Figure 4B,

Thus, ΔG = −RT ln K corresponds to the excess free energy associated with forming an equilibrium distribution of aggregates, relative to the distribution that would have been present in the corresponding random mixture. In other words, ΔG = 0 implies that the actual aggregate size distribution is essentially identical to that in the corresponding random mixture, while ΔG < 0 corresponds to aggregation and ΔG > 0 corresponds to expulsion of the ligands from the solute’s coordination shell (again relative to the corresponding random mixture). Although the above procedure for determining F



■

ACKNOWLEDGMENTS We thank Ben Widom and David Wilcox for useful comments and suggestions. We would further like to acknowledge Mike Fayer’s contribution as a role model with one foot firmly planted in the experimental domain and the other foot free to pivot strategically into theoretical realms as needed.

indicate that there is a small free-energy driving force (ΔG = −0.6 kJ/mol < 0) for forming aggregates above and beyond those in a randomly mixed solution of the same composition. The enthalpy and entropy associated with this aggregation process are characteristic of those for hydrophobic interactions, in the sense that formation of the aggregate appears to be entropicially dominated as TΔS > ΔH. Also note that the ΔH and ΔS for dissociation of the aggregates have the same (negative) sign as those for hydrophobic hydration processes at 293 K. The thermodynamic results pertaining to the aggregation of neopentane and iodide ions in water, shown in Figure 4D, indicate that aggregation is slightly disfavored (ΔG = +0.62 kJ/ mol > 0). Moreover, this aggregation process appears to be enthalpically dominated, as ΔH > TΔS. However, given the small values and larger error bars on ΔH and ΔS results, it is clear that additional studies would be required in order to reduce the statistical uncertainty of the results and, even more importantly, to determine their sensitivity to the interaction potentials and level of theory used in the simulations.

■ ■

ABBREVIATIONS MD, molecular dynamics REFERENCES

(1) Wilcox, D. S.; Rankin, B. M.; Ben-Amotz, D. Distinguishing Aggregation from Random Mixing in Aqueous t-Butyl Alcohol Solutions. Faraday Discuss. 2013, 167, online. (2) Gupta, R.; Patey, G. N. Aggregation in Dilute Aqueous tert-Butyl Alcohol Solutions: Insights from Large-Scale Simulations. J. Chem. Phys. 2012, 137, 034509. (3) Rankin, B. M.; Ben-Amotz, D. Expulsion of Ions from Hydrophobic Hydration Shells. J. Am. Chem. Soc. 2013, 135, 8818− 8821. (4) Rankin, B. M.; Hands, M. D.; Wilcox, D. S.; Fega, K. R.; Slipchenko, L. V.; Ben-Amotz, D. Interactions between Halide Anions and a Molecular Hydrophobic Interface. Faraday Discuss. 2013, 160, 255−270. (5) Pegram, L. M.; Record, M. T. Thermodynamic Origin of Hofmeister Ion Effects. J. Phys. Chem. B 2008, 112, 9428−9436. (6) Schwierz, N.; Horinek, D.; Netz, R. R. Anionic and Cationic Hofmeister Effects on Hydrophobic and Hydrophilic Surfaces. Langmuir 2013, 29, 2602−2614. (7) Rembert, K. B.; Paterová, J.; Heyda, J.; Hilty, C.; Jungwirth, P.; Cremer, P. S. The Molecular Mechanisms of Ion-Specific Effects on Proteins. J. Am. Chem. Soc. 2012, 134, 10039−10046. (8) Paterova, J.; Rembert, K.; Heyda, J.; Kurra, Y.; Okur, H.; Liu, W. R.; Hilty, C.; Cremer, P.; Jungwirth, P. Reversal of the Hofmeister Series: Specific Ion Effects on Peptides. J. Phys. Chem. B 2013, 117, 8150−8158. (9) Vazdar, M.; Pluharova, E.; Mason, P. E.; Vacha, R.; Jungwirth, P. Ions at Hydrophobic Aqueous Interfaces: Molecular Dynamics with Effective Polarization. J. Phys. Chem. Lett. 2012, 3, 2087−2091. (10) Ben-Amotz, D.; Stell, G. Analytical Implementation and Critical Tests of Fluid Thermodynamic Perturbation Theory [See Also J. Chem. Phys., 120, 4994 (2004)]. J. Chem. Phys. 2003, 119, 10777− 10788. (11) Ben-Amotz, D.; Herschbach, D. R. Estimation of Effective Diameters for Molecular Fluids. J. Phys. Chem. 1990, 94, 1038−1047. (12) Ben-Amotz, D.; Willis, K. G. Molecular Hard-Sphere Volume Increments. J. Phys. Chem. 1993, 97, 7736−42. (13) Bondi, A. Van Der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441−451. (14) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (15) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (16) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (17) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 327−341. (18) Johnson, R. D. NIST Computational Chemistry Comparison and Benchmark Database. NIST Standard Reference Database; National Institute of Standards and Technology: Gaithersburg, MD, 2011; Vol. 101

■

CONCLUSION We have presented a generalized theoretical strategy for predicting aggregate size distributions in uniformly mixed solutions and shown how such predictions may be used to both qualitatively and quantitatively analyze aggregate size distributions in molecular solutions. This strategy is illustrated using simulations of neopentane dissolved in either methane, aqueous methanol, or NaI solutions. The results clearly reveal that the significant aggregation of methane gas around neopentane at 293 K decreases markedly at 1000 K, approaching random mixing predictions. More interestingly, the present simulation results indicate that there is less aggregation of neopentane with both aqueous methanol or iodide ions in water than there is with methane in the gas phase. Methanol shows a weak tendency to hydrophobically aggregate, while iodide is weakly expelled from the first coordination shell of neopentane. We have further shown how quantitative comparisons of simulation and random mixing results may be used to obtain the partition coefficient pertaining to the adsorption of molecules from a bulk solution onto the surface of a solute and thus also the corresponding free energy, enthalpy, and entropy. The analysis of the above representative processes serves to illustrate how random mixing statistics can be used to classify and quantify other aggregations and self-assembly processes of biological and environmental interest.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

D.B.A. contributed primarily to the development of the random mixing model, and B.M.R. contributed primarily to producing and analyzing the simulation results. Both authors contributed to preparing the manuscript and have given approval to the final version of the manuscript. Funding

This work was funded by the NSF (CHE-1213338) and a PRF Research Grant from Purdue University. Notes

The authors declare no competing financial interest. G



Article

(19) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (20) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (21) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (22) Ben-Amotz, D.; Widom, B. Nonideal Gas Solvation Thermodynamics. J. Chem. Phys. 2007, 126, 104502/1−104502/10. (23) Wick, C. D.; Kuo, I. F. W.; Mundy, C. J.; Dang, L. X. The Effect of Polarizability for Understanding the Molecular Structure of Aqueous Interfaces. J. Chem. Theory Comput. 2007, 3, 2002−2010. (24) Netz, R. R.; Horinek, D. Progress in Modeling of Ion Effects at the Vapor/Water Interface. Annu. Rev. Phys. Chem. 2012, 63, 401−418. (25) Stern, A. C.; Baer, M. D.; Mundy, C. J.; Tobias, D. J. Thermodynamics of Iodide Adsorption at the Instantaneous Air− Water Interface. J. Chem. Phys. 2013, 138, 114709.

H


Analysis of Molecular Aggregation Equilibria Using ... - ACS Publications

Recommend Documents