Article Outline

Article Information

PubMed

Related Articles

  • …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 93, Issue 10, 3392-3407, 15 November 2007

doi:10.1529/biophysj.107.114181

Biophysical Theory and Modeling

Molecular Basis of the Apparent Near Ideality of Urea Solutions

Hironori Kokubo*, Jörg Rösgen†, D. Wayne Bolen† and B. Montgomery Pettitt*, Go To Corresponding Author,  

* Department of Chemistry, University of Houston, Houston, Texas
† Department of Biochemistry and Molecular Biology, University of Texas Medical Branch, Galveston, Texas

Address reprint requests to B. M. Pettitt, Tel.: 713-743-3263.

Abstract

Activity coefficients of urea solutions are calculated to explore the mechanism of its solution properties, which form the basis for its well-known use as a strong protein denaturant. We perform free energy simulations of urea solutions in different urea concentrations using two urea models (OPLS and KBFF models) to calculate and decompose the activity coefficients. For the case of urea, we clarify the concept of the ideal solution in different concentration scales and standard states and its effect on our subsequent analysis. The analytical form of activity coefficients depends on the concentration units and standard states. For both models studied, urea displays a weak concentration dependence for excess chemical potential. However, for the OPLS force-field model, this results from contributions that are independent of concentration to the van der Waals and electrostatic components whereas for the KBFF model those components are nontrivial but oppose each other. The strong ideality of urea solutions in some concentration scales (incidentally implying a lack of water perturbation) is discussed in terms of recent data and ideas on the mechanism of urea denaturation of proteins.

Introduction

Molecules in cells function in a highly crowded, concentrated, nonideal solution environment. Therefore, the usual treatments that make use of concepts from the ideal, infinite dilution solution limit are quantitatively inadequate for many biological problems. Urea is an interesting case. It has a strong concentration-dependent effect on the folding/denaturing transition for proteins in solution. This might appear to be due to nonideality in solution. In fact, in some concentration scales it shows substantial deviations from ideality. Yet, in the molar scale, it appears almost ideal. Much literature has been devoted to how urea must change water's structure to make its solution such a powerful denaturant 1. Theoretical work from this lab and experiments from others question the water structure change hypothesis 2,3.

The effect of osmolytes and a multicomponent environment in general is to change the chemical potentials of the components, most notably that of the macromolecular solutes. Such changes result in differences in stability for conformations and oligomerization versus simple aqueous solutions. The effects can be quite significant for biological systems where common use is made of urea, proline, sucrose, glycerol, etc. 4,5. To consider these effects we must consider the system's free energy (or activity) in general.

The Gibbs free energy change, dG, is

(1)
where S is the entropy, T is absolute temperature, P is pressure, ni is the number of the ith species and μi its chemical potential. The systems in which we are interested, for the most part, in this work are isothermal-isobaric systems. In a multicomponent mixture, the thermodynamic condition of the system is expected to be strongly affected by the chemical potentials. However, interpretation and even measurement of chemical potential or activity is confounded by standard state and even the concentration scale used.

It is convenient to start from the understanding of ideal solution properties to understand real systems. Ideal solutions are a convenient if occasionally misleading approximation to real solutions, much as an ideal gas is an approximation to real gases. There are two well-known ideal solutions. One is the symmetric ideal solution and the other is a dilute ideal solution 6. Isotopic and some isomeric mixtures are typical examples which are often very close to a symmetric ideal solution. However, most pure electrolytes and osmolytes are in the solid state for the range of temperature and pressure in which we are interested. It is notoriously difficult to determine the relative activity of the solid state with respect to the pure liquid experimentally. In such cases, the infinitely dilute state can often be used as the standard state, and one then considers the deviation from the dilute ideal solution.

No real solution is an ideal solution exactly. Adding solutes to liquid water generally causes, to a greater or lesser extent, colligative effects such as vapor pressure lowering, boiling point elevation, freezing point depression, and osmotic pressure changes. These properties depend only on the number of solute molecules in the case of a dilute ideal solution. They are relatively independent of chemical properties of solutes for dilute solutions. They generally depend on species and concentrations of solutes in real solutions.

The existing solution theories which use the infinitely dilute state as the standard state most often explain solution behavior only in a very dilute solution. There is no complete, analytic theory which can explain experimental data at concentrations of >1.0mol/L. In development of such theory, Rösgen et al. were recently very successful in fitting partition function ratios as parameters to experimental data in an activity series with analysis based on a derivation arising from a semi-grand canonical ensemble 7,8. Here, using urea solutions as a case of interest, we would like to address whether the fitted parameters are true representations of the ratios of low order partition functions or just effective coefficients, i.e., direct representations or renormalized. We discuss the dependence on the urea model parameters as well in Results and Discussion.

It is generally cumbersome to obtain experimental solvation free energy (excess chemical potential) in most electrolytes and osmolytes. One reason is that relative activity coefficients depend nontrivially on the concentration scale used. Another reason is because the vapor pressure at room temperature is often too small to measure accurately. Smith et al. 9 estimated the solvation free energy of urea computationally. However, their estimation is based on some assumptions. Because there is no experimental data on the vapor pressure of urea at room temperature (urea is essentially nonvolatile), they estimated from an extrapolated value of the vapor pressure at high temperatures. It is difficult to validate these estimations because existing experimental data on urea vapor pressure is inconsistent, and with questionable accuracy 10,11. The validity of the extrapolation method to low temperatures is similarly unclear. Urea solutions show a striking apparent dependence on how we view them in terms of standard state and concentration scale. They are close to a dilute ideal solution in the molarity scale, but not for other scales such as the mole fraction or molality, nor for a symmetric ideal solution 9 as we show below. We also show that it is possible to understand the near ideality of urea solutions that occurs with certain scales and reference systems by molecular simulation with appropriate theoretical analysis.

There is no direct experimental solvation free energy (excess chemical potential) data for urea solutions. Thus, it is difficult a priori to judge the applicability of the various force fields for urea from the solvation free energy value calculated using simulations. Experimental aqueous solvation free energy data exists only for a few solutes such as some small alkanes, alcohols, and amides, because such solutes are volatile at room temperature. Instead, for urea we have only the experimental activity coefficient data 12,13,14. Experimental activity coefficients are obtained from osmotic coefficient data on urea with knowledge of the osmotic pressure of water via the Gibbs-Duhem relation.

In principle, for any force field, the activity coefficient of urea can be obtained through calculation of the chemical potential at different concentrations, although that requires high precision to compare with experiment. In this article, we examine how and why the activity coefficient changes in different concentrations of urea solutions. We use free energy calculations with sampling from molecular dynamics simulations with two different all-atom force fields to generate hypotheses about the mechanism to explain the activity at the molecular level. We hope to clarify the origin of the behavior of urea in aqueous solution as a prelude to considering its profound effect on proteins.

Calculating the chemical potential with sufficient precision to obtain activity coefficient changes is still a challenging computational problem. It is necessary to calculate the chemical potential quite precisely to estimate the change in an activity coefficient by simulations. Here we combine and contrast the thermodynamic integration method, perturbation method, and Bennett's acceptance ratio method for our activity coefficient calculations. We compare simulation results with experimental data over a wide range of urea concentrations.

In the next section, we review the theoretical framework of ideal and nonideal solutions as well as the often confused concentration dependence. The calculational methods to obtain sufficiently precise chemical potentials are also explained. Readers familiar with the connections of Raoult's and Henry's laws to the modern theories of solution may skip Ideal Solutions. In Methods, the models and details of the simulation are presented. In Results and Discussion, in the context of both the models, the data are given. Conclusions is devoted to our remarks about what our results imply about the mechanism of action of urea on proteins.


Theory

For our subsequent analysis, we must separate and quantify the effects of standard state from concentration-scale dependence of the chemical potential or activity coefficients. We start with a discussion of ideality which has more than one definition in the literature. This will be important to relate molecular level correlations to the various measures of nonideality.

Ideal solutions

We first review the theoretical definitions of an ideal solution to set our work in context. The reader familiar with this area and interested in the case of urea solutions may skip to Methods. The concept of an ideal solution was first developed by Raoult historically 15. Raoult found that the vapor pressure of a component over a liquid solution at equilibrium is proportional to the mole fraction of the component in solution,

(2)
(3)
where is chemical potential of the substance A in the mixture solution, is chemical potential of A in pure A, pA is the partial vapor pressure of A in the mixed solution, and is the vapor pressure of pure A.

Here, the ideal solution is defined as the solution which satisfies the following relation for mole fractions 0≤xA≤1:

(4)
In this definition, it is not necessary to assume that vapor is an ideal gas. This type of ideal solution has been called a symmetric ideal solution 6. In Eq. (4) the second term, right side, is always ≤0.0, and it implies that the chemical potential in a mixture is always smaller than that in the pure liquid. Note that symmetric ideal solutions are defined at constant temperature and pressure.

The Gibb's free energy difference per molecule or chemical potential, Δμ, when an ideal solution is made from nA of substance A and nB of substance B becomes

(5)
with change in enthalpy, ΔH, and entropy, ΔS,
(6)
(7)
Thus, there is no enthalpy change (no heat of mixing) (Eq. (6)). Equation (7) is the entropy of mixing in an ideal gas. Therefore the difference in chemical potential between the system of nA moles of pure A and the mixture system of solution of nA moles of A and nB moles of B is
(8)
This makes clear that the second term on the right-hand side of Eq. (4) is the contribution from the entropy of mixing.

What are necessary and sufficient conditions for such ideal solutions? By differentiating Eq. (4) by xA we have

(9)

On the other hand, we know from Kirkwood-Buff theory 16,

(10)
where ρ is the number density, ρ=ρA+ρB=〈NA〉 /V+〈NB〉 /V, and
(11)
is called the Kirkwood-Buff integral which requires the radial pair distribution, gαβ(r), as a function of distance, r, between molecular species α and β.

Therefore, the necessary and sufficient condition is that there are no excesses or deficits of one molecule around another versus a simple random distribution on average or that

(12)

In many solutions, the vapor pressure of solvent of very dilute solute solutions obeys

(13)
where, from now on, we consider A as solute and B as solvent. This is the empirical law for infinitely dilute solutions at constant T. Thus, in particular, pure solvent is an ideal solution. When the solvent obeys Raoult's law in dilute solutions, the solute follows
(14)
Here, pA is vapor pressure of solute and kA is a constant which depends on the species of the solute. This relation is derived straightforwardly by integrating the Gibbs-Duhem relation. This is called Henry's law 17. In this case, the chemical potential of solute becomes
(15)
at xA∼0 (we assume that the gas phase is an ideal gas or very dilute). Henry's law can be derived differently as follows. When the mole fraction of solute approaches 0, its vapor pressure also approaches 0, and so
(16)
is another statement of Henry's law. Therefore, the existence of this proportionality relation does not depend on the species of the solutes or the number of components of the solution. The proportionality constant kA contains the actual dependence on the substance. Similarly, we obtain the following relations by using different units such as molality, m, or molarity, ρ, instead of mole fraction, respectively,
(17)
(18)

The vapor pressure of the solute in an infinitely dilute solution is proportional to mole fraction, molality, or number density. Here, we showed the proportional relations to vapor pressure, but one can derive similar relations for activity.

The validity of the concentration dependence of Henry's law depends on the species of the substances involved, but is often a good approximation at low concentration where we do not expect significant associations or influences among the solute molecules. If Henry's law is satisfied at high concentration, the chemical potential is the same as that of the infinitely dilute solution, that is, a dilute ideal solution.

If the chemical difference between solute and solvent is very small and the solution is a symmetric ideal solution, Henry's law is then satisfied as long as we use mole fraction or molarity scales. In other words, if it is a symmetric ideal solution with similar molecular volumes, it is also a dilute ideal solution in mole fraction scale and molarity scale. The dilute ideal solution in the molality scale is exceptional and it is not a symmetric ideal solution even if it is a dilute ideal solution as explained later.

The regular solution defined by Hildebrand is the solution which satisfies the following 18:

(19)

This differs subtly from Raoult's definitions in Eqs. (6) and (7) by the allowed change in ΔH. This is consistent with the assumption that there are no specific interactions such as associations between molecules so that the distribution and orientation of molecules are completely random. Therefore, a urea solution is not expected to be regular because it has specific strong interactions and local orientational preferences via hydrogen bonds. In regular solutions, thermal fluctuations are assumed to be strong enough to overcome the specific or orienting interactions between different molecules and cause random mixing. If the difference between two molecules is large and random mixing does not occur, solubility is generally small in this case and it becomes a very dilute solution. As a result, there will be little direct interactions between solute molecules.


Concentration scale dependence of nonideality

We now discuss general solutions which deviate from ideal solutions and the dependence of measured changes in activity on concentration scales. We may write

(20)
Here, aA is relative activity and γA is the activity coefficient. The value is the chemical potential of pure A, so γA=1 when xA=1. Various nonideal effects are all included in γA in this expression. Note that aA is deemed “relative” because we used a standard state. If we consider μA without choosing a standard state,
(21)
where fA is the absolute activity.

As above, a natural choice of the standard state is a pure solute solution. However, most pure salts and many osmolytes are solid and not liquid in the range of temperature and pressure in which we are usually interested (∼1atm, ∼298K). In this case, which is the chemical potential of the pure liquid, is not measurable experimentally. Another natural choice of the standard state is the infinitely dilute solution. Therefore, in many cases we take infinite dilution as the standard state for osmolytes.

The fact that relative activity is thus defined in various, disparate ways causes confusion in the literature. Only when the relative activity with respect to the pure liquid can be evaluated experimentally is it the custom to take pure liquid as the standard state. We often take the infinite dilute solution as the standard state because it can be referred to easily experimentally.

The activity coefficient at mole fraction x when the standard state is infinitely dilute may be introduced in the following way:

(22)
Here, is the activity coefficient of the infinite dilute solution. Using Eq. (20) becomes
(23)
using the infinite dilution standard state. Thus, we define the sum of the chemical potential of pure A and the contribution from the infinitely dilute standard state as
(24)
which can be referenced experimentally.

We next consider the consequences of concentration scale change from mole fraction to molality [mol/kg]. Molality can be expressed

(25)
where MB is the mass [g/mol] of the solvent B, and nA and nB are the mole numbers of the substances A and B in the system. The relation between mole fraction and molality is simply
(26)
Substituting Eq. (26) in the first equation of Eq. (23), we have
(27)
where m0=1 [mol/kg] was introduced to make the content of the ln dimensionless and
(28)
With this we define
(29)
Equations (28) are the activity coefficient and chemical potential in the infinite dilution standard state with the molality scale, respectively.

Similarly for molar concentration units referenced to infinite dilution,

(30)
where [mol/L] is the number density of the pure solvent B and again makes the argument of the dimensionless and makes the activity coefficient in the infinitely dilute solution 1.0 (see Eq. (31) below). The activity coefficient and standard state chemical potential in the molarity scale are
(31)
(32)
(33)

The well-known relation between the molality and the molarity activity coefficients for a given temperature and density can be obtained similarly:

(34)

Alternatively, one could just use the analytic transformation between the concentration scales, which does not require knowledge of the density at the relevant temperature previously derived, to good accuracy 7.

The chemical potential of solute A in ideal dilute solutions (xA ∼ 0, xB ∼ 1) are thus expressed in any of the following ways,

(35)
and the chemical potentials in terms of activities, and activity coefficients for general nonideal, dilute solutions are
(36)
The molality and mole fraction of solutions do not depend on temperature, but molarity does depend on temperature because the volume and thus, density, changes with respect to temperature. Given this and the volume-versus-mass issue, the apparent nonideality of a given solution is qualitatively very different if we choose different concentration scales except in the limit of infinitely dilute solutions where γ becomes 1.0 in every scale.


Reference state forms of chemical potential

In this section, we formally explore the consequences of interpreting the chemical potential changes for different reference systems in various concentration units. The object here is to phrase the relevant relations in terms of quantities readily computable from simulation or liquid state theory. Besides the numbers of molecules of each species, we will require the average volume, 〈V〉 and the excess chemical potential.

We perform our simulations found in Results and Discussion in the isobaric-isothermal or NPT ensemble for the calculations of the excess chemical potential. Equation (59) could be used to calculate chemical potential but refers to a different ensemble. In general, one could simply Legendre-transform the results. However, because the correlation between volume and energy in Eq. (59) is small in urea solutions (the correlation coefficient magnitudes in our simulations were ∼0.15), we can reliably, approximately transform Eq. (59) to

(37)
Here, the ensemble average is taken over and
(38)
Urea solutions at room temperature and pressure are dense liquids and urea is too large for a successful implementation of the particle insertion method for We evaluated three other methods below, thermodynamic integration, perturbation theory, and the Bennett-Pande acceptance ratio method, for an estimation of of Eq. (37). This term is often interchangeably referred to as either the excess chemical potential or the solvation free energy. The detailed conditions of our simulations are described later in Methods.

Given the breakdown of μA in Eq. (37), we now write the chemical potential form in terms of quantities readily available from simulation in the mole fraction, molality scale, and molarity scale below:

(39)
where the terms with * are the values at the infinitely dilute state.

Let us first consider infinite dilution as the standard state. Consider the deviation from the dilute ideal solution. In the mole-fraction scale standard state, the chemical potential and activity coefficient become

(40)
(41)

In the molality scale, we thus have

(42)
(43)

We see that the relationship between activity coefficients in the mole fraction and molality case is the same as shown in Eq. (28). Namely, Eqs. (41) and (43) are consistent with Eq. (28) as they should be.

In the case of the molarity scale, however, we have

(44)
(45)

These equations (Eqs. (41), (43), and (45)), in terms of readily computable quantities, show how the meaning of nonideal solutions change when the concentration scale changes. The activity coefficient in the molarity scale measures only the excess free energy difference (solvation free energy difference). As is well known but often underappreciated, the apparent deviation from ideality for dilute solutions strongly depends on the scale.

We next consider the nonideal deviation for a symmetric ideal solution reference. Pure solute is the standard state in this case. However, symmetric ideal solutions are usually defined with the mole fraction scale (Eq. (4)), not with molality or molarity. Therefore, the standard state and activity coefficient may be written as

(46)
(47)

where the terms with ** are the values for a pure solute state. There are two ways to make the activity coefficient 1.0:

Case 1:
Case 2:
In a mixture of very similar molecules, case 1 is approximately achieved. This is the familiar example of a symmetric ideal solution. Case 1 also makes Δ in Eq. (41) and Δ in Eq. (45) 0.0. Therefore, if it is a symmetric ideal solution and the number of total molecules in the system per volume is constant (), it is also a dilute ideal solution in mole fraction scale and molarity scale. However, in Eq. (43) is not 0.0 because the number of solvent molecules changes, so symmetric ideal solutions are not dilute ideal solutions in the molality scale.

We may also take the pure solute as the standard state in molality scale and molarity scale. This choice, however, does not measure the deviation from symmetric ideal solutions because symmetric ideal solutions are usually defined by the mole fraction scale. In the molality scale, the standard state and activity coefficient may be defined as

(48)
(49)
Without any solvent molecules, we have a numerical problem. In Eqs. (48) and (49) the standard state in pure solute effectively causes
(50)
Therefore, it is not reasonable or convenient to take the pure solute state as the standard state in the molality scale.

For the molarity scale, the standard state and activity coefficient would be

(51)
(52)
Clearly, standard state and concentration scale choices change the meaning of the activity coefficient in a qualitative and quantifiable way.


Excess chemical potential calculation by simulation

We next require an accurate way to calculate the change in excess free energies or chemical potentials in solution for adding a solute. To explore mechanism and to control for force field, two variants are used below, well-known OPLS 19 and the newer KBFF by Smith and co-workers 20. We demonstrate the efficiency and precision characteristics for three methods of calculation. We present a brief review of the methods here first for coherence. Considerably more detailed technical reviews exist in the recent literature 21.

Thermodynamic integration method

The well-known thermodynamic integration method calculates the free energy difference between the state i and the state j by

(53)
Here, we assumed the potential energy function U(λ) is written as a function of a coupling parameter, λ, and λ=0 and λ=1 correspond to the state i and the state j, respectively. We can define many different functional dependencies on λ corresponding to different integration paths. The simple linear ramp is
(54)
We consider the beginning and final system to correspond to N particle system and N+1 particle system, respectively. There is a well-known pole at the origin for the Lennard-Jones interaction with respect to λ. To avoid numerical instability, a nonlinear function is sometimes used to alleviate large absolute values at a small λ. The following soft-core potential function was developed to avoid singularity in van der Waals (vdW) interactions, to calculate the free energy difference precisely 22,23:
(55)
We adopted this soft-core potential in our calculations.


Perturbation method

In so-called thermodynamic perturbation methods we derive without any approximations:

(56)
This equation means that we can obtain the free energy differences between the state i and the state j by calculating the ensemble average of exponential of the potential energy difference at the state i ensemble. This method is only useful when the state j may be conveniently sampled from state i. To avoid difficulty one typically divides the region into many subregions:
(57)
(58)
This method of small windows or steps is exacerbated if the energy barrier in or between subregions is large. A different approach to a solution is to use nonphysical sampling such as embedding the problem in a higher dimensional space or using generalized ensemble methods. We found we could use Eq. (57) or (58) because our urea solution was not so complicated.


Widom test particle insertion method

For completeness we mention the Widom insertion method which has uses both conceptual and practical 21,24. The basic equation of Widom method 24 can be derived like Eq. (56) in the NPT ensemble,

(59)
where we show the equation in the case of a monoatomic. In the NVT ensemble we have
(60)
In this method we insert a particle randomly in the system and calculate the potential energy with which the inserted particle interacts. This method works very well in low density systems, but fails in dense liquids and solids especially when the inserted molecule is large. Therefore the precision of calculation is often not sufficient, especially in calculating activity coefficients, with this method.

In the case of polyatomic molecules, the chemical potential becomes the following one instead of Eq. (59),

(61)
where q represents the rotational, vibrational, electronic, and nuclear partition function terms. We assume that these terms can be separated from the configuration integral without concern about whether they are expressed classically or quantum-mechanically. This assumption breaks down when strong intramolecular interaction changes the counting of degree of freedom. In this article we set q=1 (Eq. (59)).


Bennett-Pande ratio method

The Bennett method 25 for optimizing sampling was originally implemented to accelerate convergence of Monte Carlo methods. More recently Pande and co-workers 26 used this principle to achieve optimal sampling in a variant of the familiar thermodynamic perturbation theory method. The input data to calculate the free energy is the same essentially. The difference with perturbation methods is the numerical precision. Bennett's method minimizes statistical error. The two basic equations are

(62)
where
(63)
Here, n0 and n1 are the sample numbers in the ensembles at the state 0 and the state 1, respectively. In practice, we plot both sides of Eq. (62) as a function of C, and solve for the C which satisfies this equation. We then have
(64)

If we collect the same number of samples (n0=n1), we can calculate the free energy difference:

(65)
This method uses the same input data as perturbation method (U1–U0), but recent studies show that it is often the best way to obtain free energy differences 26,27,28.




Methods

We wish to understand the mechanism by which a model might reproduce experiment. As a control and illustration, we evaluated two different urea models, namely OPLS 19 and KBFF 20 and examine the dependence of the activity coefficients on the force-field parameters. We show the parameters of these two models in Table 1. The geometry of the urea models is the same. We adopted TIP3P model for water and the minor consequences have been noted previously 2.

Table 1 Force-field parameters for the OPLS model 19 and KBFF model 20
MassChargeɛ(kJ/mol)σ(nm)
OPLS
O15.999−0.3900.878640.296
C12.0110.1420.439320.375
N114.007−0.5420.711280.325
H111.0080.3330.000000.000
H121.0080.3330.000000.000
N214.007−0.5420.711280.325
H211.0080.3330.000000.000
H221.0080.3330.000000.000
KBFF
O15.999−0.6750.560000.310
C12.0110.9210.417000.377
N114.007−0.6930.500000.311
H111.0080.2850.088000.158
H121.0080.2850.088000.158
N214.007−0.6930.500000.311
H211.0080.2850.088000.158
H221.0080.2850.088000.158

KBFF urea model was developed to reproduce the experimental activity data. In Weerasinghe and Smith 29, they determined a charge distribution for urea atoms by using Kirkwood-Buff integrals obtained from simulations. Kirkwood-Buff relations yield the derivative forms of activity coefficients. It is necessary to integrate them subsequently by, for instance, assuming the experimentally suggested functional form. However, Kirkwood-Buff G factors converge very slowly and it is difficult to obtain precise values from simulation.

The OPLS parameters have been used extensively in the literature to model and simulate a variety of systems. The OPLS model was fit to several liquid state properties including heats of solvation among others 19.

By using two different models we hope to get an idea of the force-field dependence in the implied mechanism. In this article, for these potentials, we calculate the chemical potential directly from the simulations by the methods described in Theory and thereby obtain the relative activity coefficients to compare directly with experiment. Thus, we test multiple methods over a range of concentrations in hopes of obtaining sufficient precision to evaluate the accuracy of the models for this purpose. In future work we will consider three component solutions which include a polypeptide.

We prepared urea aqueous solutions at seven different concentrations for the OPLS and KBFF models from dilute solution to the pure urea. We first estimated the molecular volume for one urea roughly to obtain the expected concentrations. It was approximately two and half times the volume of one water. For the pure solute solution state, since urea is a solid at room temperature, a pure urea sample was equilibrated at a higher temperature and super-cooled to 298K.

To prepare our solutions we took an equilibrated water box and the urea box, and randomly removed the required urea and water molecules from each of these boxes in turn to achieve the required numbers for the desired solution concentrations. Rather than the standard replacement, we put these two boxes in contact with the normal periodic boundary conditions at the large volume equivalent to the two pure liquids. Next the volumes were shrunk to the target value to mix the urea and water. Thus, we prepared the systems by deciding the number of urea and water molecules and shrinking to a given volume. We performed a minimization for 500 steps and mixed for 100ps of NVT simulation for each concentration. Temperature was controlled by Nosé method to 298K 30. In this process it was found that the mixing occurred not only spontaneously but in fact quite rapidly, on the order of the time to shrink the box. We next performed 300ps NPT equilibrations at 298K and 1atm using the Nosé-Anderson method for the temperature and pressure control 30,31. The final system sizes are close to 34Å×34Å×34Å in each case.

We thus obtained the initial configurations of urea solutions at seven different concentrations for two different urea models. The number of total systems considered is 14 (7 concentrations×2 urea model). All the systems in our simulations are listed in Table 2. We show typical snapshots of configurations of KBFF urea solutions in Fig. 1, which qualitatively confirms the good mixing.

Table 2 Symbols are Nu, number of urea molecules, and Nw, number of water molecules
NuNwMole fractionMolalityMolarityVolumeDensity
OPLS113050.00076570.042530.0428838.731.011
4711880.038062.1962.02938.471.046
9510770.081064.8964.08938.581.081
1429550.12948.2546.15238.321.115
1898380.184012.528.20938.231.149
2486610.272820.8311.0437.311.193
53001.0—22.6538.851.361
KBFF113050.00076570.042530.0428638.751.010
4711880.038062.1962.01338.781.037
9510770.081064.8964.02739.171.064
1429550.12948.2546.03139.101.093
1898380.184012.527.99339.271.118
2486610.272820.8310.7038.501.156
53001.0—22.4139.271.346
The units of molality, molarity, volume, and density are (mol/kg), (mol/L), (nm3), and (g/cm3) respectively.

Because all systems have very similar box size, we used the same cutoff length 15Å for vdW interactions, which was >4σ. We used a link-cell Ewald method for electrostatic interactions 32. In the studies undertaken here it may be quite important to estimate the electrostatic interactions accurately to correctly distinguish among the potential energy models. Free energy calculation is known to be very sensitive to the method used for electrostatic energy calculation 33.

We calculated the chemical potentials by Eq. (37) inserting one urea molecule. The term for of Eq. (38) was calculated using the temperature and average volumes. of Eq. (38) was calculated by three different methods. Those are the thermodynamic integration method, perturbation method, and Bennett acceptance ratio method. In calculating we divided the potential energy of the inserted urea molecule into the vdW and electrostatic terms for subsequent analysis.

In the case of vdW interactions, a soft-core potential (Eq. (55)) was used. The range for λ was divided into 50 subregions (= 51 points). We thus calculated 51 λ-points of the integrand of Eq. (53) for the thermodynamic integration method and for the perturbation method using Eqs. (57) and (58). The same sampling data as used for the perturbation method was used for the Bennett method after Pande. In the case of the electrostatic interaction contributions, we used Eq. (54) and divided λ into 25 subregions (= 26 points). By performing the repulsive core first, there is no remaining singularity, so it is not necessary to use soft-core potential for electrostatic interactions.

In the case of thermodynamic integration method, plotting the integrand of Eq. (53) at the calculated λ-points and integrating numerically we obtain the excess chemical potential. In the case of the perturbation and Bennett methods, we used Eqs. (57), (58), and (62), respectively. Summing up the subdivisional free energy difference at the calculated λ-points produces the excess chemical potential. The input data to Eq. (62) is Uj–Ui and this is the same as the perturbation method case. As we show in the section below, these three different methods give us almost the same values if the sampling yields enough precision as expected. We show some example figures of the free energy calculation paths below.


Results and discussion

Chemical potential calculations

We now consider the results of the calculations of μexcess as well as μ as a function of concentration and standard state. Figure 2a shows the integration path in the case of the most dilute KBFF solution. The error bar of each point was estimated by dividing the data into 10 blocks and calculating the standard deviation of the average values. The error bars were largest in the region of most curvature. That region demands more sampling for precise estimations. The points between 0.26<λ<0.60 in Figure 2a were calculated from 1-ns simulations, and other points were from 160-ps simulations. Integrating this path from 0 to 1 gives the free energy difference of the vdW transfer from vacuum to solution. Figure 2c is the corresponding figure for the electrostatic part. Summing up the integrated values of Figure 2ac, becomes the total free energy difference (the excess chemical potential).

Display large version of this figure
Figure 2
The excess chemical potential integrand components in the case of the most dilute KBFF urea solution. (a) Integration path for the calculation of vdW terms of excess chemical potential by thermodynamic integration method. Integrating this path about lambda becomes the total vdW excess chemical potential. (b) The 51 λ-points for the calculation of vdW terms of excess chemical potential by Bennett acceptance ratio method (black asterisk) and perturbation method (+, estimation Eq. (57), blue circle; –, estimation Eq. (58), red triangle). Adding the differences of 51 points becomes the total vdW excess chemical potential. Asterisk and circle overlap very well. (c) Similar to panel a for electrostatic terms. (d) Similar to panel b for electrostatic terms.

Figure 2bd, show the calculated values of the λ-points at the most dilute solution for the perturbation method and Bennett method using the KBFF model. We see that Bennett estimation points almost overlap perturbation ones. Summing the values yields the total free energy difference.

Based on the experimental trends we expected the deviations from ideality to be difficult to obtain with sufficient precision to evaluate the difference between models. Thus we tested the precision and convergence of the three different free energy techniques. We confirmed that the free energy values calculated by these three different methods gave similar values for each system. For example, in the system of Fig. 2 the obtained values by thermodynamic integration method, perturbation method, and Bennett method are respectively 2.147, 2.113, and 2.138kJ/mol for the vdW part, and −46.411, −46.393, and −46.398kJ/mol for the electrostatic part, respectively. Testing of convergence in model systems showed that the Bennett method as implemented by Pande and co-workers 26 gave the least errors for a given amount of sampling (data not shown). Given that and the subkJ/mol precision of all the methods, we used the estimations by the Bennett method in the following analysis for all the other systems.

Table 3 shows the chemical potentials and their components obtained from our simulations. Because μA of OPLS is more negative than that of KBFF at the same concentration, OPLS urea dissolves in aqueous model solutions using the TIP3P water model better than KBFF urea.

Table 3 Chemical potential and its components for OPLS and KBFF model urea solutions using Eq. (37)
Mole fractionμA
OPLS
0.0007657−41.310.87−56.00−55.13−96.44
0.03806−31.760.66−56.40−55.74−87.49
0.08106−30.020.33−56.19−55.86−85.88
0.1294−29.010.49−56.56−56.07−85.08
0.1840−28.290.040−56.63−56.59−84.88
0.2728−27.56−0.63−57.11−57.73−85.29
1.00−25.78−1.29−56.43−57.72−83.50
KBFF
0.0007657−41.312.14−46.40−44.26−85.57
0.03806−31.781.29−45.57−44.29−76.06
0.08106−30.060.38−44.79−44.41−74.46
0.1294−29.06−0.27−44.22−44.49−73.54
0.1840−28.36−1.16−43.12−44.28−72.64
0.2728−27.64−2.44−41.83−44.28−71.91
1.00−25.80−4.82−37.40−43.23−69.03
The quantity consists of two terms, the vdW part and the electrostatic part plus becomes in Eq. (37). The quantity plus becomes the chemical potential μA. The units of chemical potential are (kJ/mol).

The ideal part of the chemical potential does not depend as strongly on the interaction model (see, for instance, Eq. (38)). In fact only small differences in volume or density are model-dependent. The deBroglie wavelength term is clearly common for OPLS and KBFF solutions. increases as the concentration increases because the number of urea molecules per volume increases (see Eqs. (37) and (38)). Thus, the entropy at higher urea concentration (mole fraction) is smaller than in the lower one, as it should be.

We see that of KBFF solutions, the excess chemical potential, is almost constant except for the pure urea system. This requires a remarkable cancellation to obtain the same excess solvation free energy at different urea concentration for the KBFF force field. On the other hand, for the OPLS force field, decreases as the concentration increases. For such a force field, we see that the total chemical potential change and the excess part move in different directions. Interpreting only the excess solvation part of the free energy, as is often done in simple modeling, would indicate that urea dissolves in higher urea concentration solutions more favorably.

In Table 3 we examine the potential components of that is, and </