3052
Ind. Eng. Chem. Res. 1998, 37, 3052-3057
A Theoretically Improved Perturbation Model for Activity Coefficients of Amino Acids and Peptides in Aqueous Solutions Mohammad K. Khoshkbarchi Hyprotech Ltd., 300, 1110 Centre Street North, Calgary, Alberta, Canada T2E 2R2
Juan H. Vera* Department of Chemical Engineering, McGill University, Montreal, Quebec, Canada H3A2B2
A theoretically improved model, based on the perturbation theory, to correlate the activity coefficients of amino acids and peptides in aqueous solutions is presented. The radial distribution function for the reference system and for the perturbation term have been obtained from the solution of the Percus-Yevick equation for hard sphere systems. The pair interaction energies required by the model are those due to dispersion forces and dipole-dipole interactions. These interactions are represented by a Lennard-Jones (6-12) and a Keesom expression, respectively. The depth of potential well in the Lennard-Jones model and the size parameter of hard spheres were evaluated by the minimization of the Helmholtz free energy of the systems. The significance of the parameters is justified by comparison of the model with the Gibbs-Bogoliubov inequality. Independently evaluated dipole moments of the amino acids are used. The model can accurately correlate the activity coefficients of amino acids and peptides in aqueous solutions over a wide range of concentration. Introduction
Thermodynamic Framework
Models for the activity coefficients of amino acids, based either on empirical local composition or group contribution methods, have been recently discussed on a comparative basis by Khoshkbarchi and Vera (1996a; 1996b). On a line closer to this work, we recently presented a model based on the perturbation theory to correlate the activity coefficients of amino acids and peptides in aqueous solutions (Khoshkbarchi and Vera, 1996b). This model had the advantage of a simple formulation obtained using the simplifying assumption of a uniform radial distribution function in the perturbation term. In the present study we have improved the theoretical structure of the perturbed hard sphere model considering for both the reference system and the perturbation term the radial distribution function obtained from the solution of the Percus-Yevick (P-Y) equation for hard sphere systems (Lebowitz, 1964; Blum and Høye, 1977). Due to the large dipole moment of amino acids (Cohn and Edsall, 1965), the interaction energies for the perturbation term are considered to be due to dispersion forces, represented by a Lennard-Jones (6-12) expression (L-J), and to dipole-dipole (D-D) interactions, represented by an angle-averaged D-D interaction in the form of the Keesom equation (Maitland et al., 1981). The depth of potential well and the size parameter of amino acids and peptides were evaluated by minimizing the difference between the activity coefficients calculated from the model and experimental data obtained from the literature. The physical significance of the parameters is justified by comparing the model with the Gibbs-Bogoliubov inequality (Hansen and McDonald, 1976).
The perturbation theory developed by Zwanzig (1954) was applied to the canonical ensemble partition function by Rowlinson (1964) and by McQuarrie and Katz (1966). The theory was extended to liquids by Henderson (1974). In the perturbation theory, the total interaction energy of the system is divided into a contribution of a reference system and a perturbation contribution. The perturbation term is expressed in the form of a series in which the first-order term involves the two-body interactions and the higher-order terms involve the three- or higher-body interactions. This theory has been shown to be able to accurately represent the Monte Carlo simulation data for many systems, such as mixtures of L-J potential molecules (Grundke et al., 1973) and mixtures of square-well potential molecules (Bokis and Donohue, 1995). Compared with phenomenological local composition or group contribution methods, perturbation theory has the advantage of ease of incorporation of different forms of interactions arising from different structure and energy characteristics of the components in the system. To model the activity coefficient of amino acids in aqueous solutions, the perturbation theory in a primitive framework has been employed. The primitivity of the model implies that the solvent in the system is treated as a dielectric continuum and the solute/water system is reduced to the treatment of a single component system. For simplicity, the model is restricted to a single solute in an aqueous solution. An extension of the model for multicomponent systems is presented in the Appendix. The residual Helmholtz free energy of the system, Ar, is considered to be composed of a contribution from a reference system, and a perturbation term:
* To whom correspondence should be addressed.
Ar ) ARef + APer
S0888-5885(97)00806-3 CCC: $15.00 © 1998 American Chemical Society Published on Web 04/29/1998
(1)
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3053
where the superscripts Ref and Per denote the contributions of the reference system and perturbation term, respectively. In this work, the contribution of the reference system to the residual Helmholtz free energy is represented by an equation of state for hard spheres obtained from the solution of the P-Y equation for hard sphere systems (Lebowitz, 1964). This equation for a single-component system, is written as
3η(2 - η) ARef Ar,hs ) ) -ln(1 - η) + kT kT 2(1 - η)2
(2)
∫σ u(r) ghs(r) r2 dr ∞
(3)
(4)
where N is the number of the molecules in the system, u(r) is the interaction energy as a function of intermolecular distance, r, and ghs(r) represents the radial distribution function of the reference term (i.e., of the hard sphere system). To facilitate the further derivations of equations, eq 4 is written in a dimensionless form as
APer ) 2πFσ3 NkT
∫1∞u*(r*) ghs(r*) r*2 dr*
uL-J(r) 4 -12 ) (r* - r*-6) kT kT
u*(r*) ) u*L-J(r*) + u*D-D(r*)
(9)
∫1∞u*L-J(r*) ghs(r*) r*2 dr* + ∫1∞u*D-D(r*) ghs(r*) r*2 dr*]
(10)
According to the perturbation theory, the radial distribution function used in the perturbation term should be the same as the one used in the hard sphere term. However, an explicit expression for the radial distribution function of the hard sphere systems, ghs(r*), is usually not readily available and the only available information is the Laplace transform of r*ghs(r*). For the P-Y hard sphere equation of state the expression of the Laplace transform of r*ghs(r*) for a singlecomponent system reads (Lebowitz, 1964):
G(s) ) L [r*g(r*)] )
∫0∞r*g(r*) e-sr
*
dr* ) L(s) e-s
(1 - η)2Q(s) s2
(11)
where L is the Laplace transform operator and the expressions for L(s) and Q(s) are given as follows:
Q(s) )
P(s) + 12ηL(s) e-s
(12)
(1 - η)2s3
(
L(s) ) 1 +
η s + 1 + 2η 2
)
(13)
P(s) ) (1 - η)2s3 + 6η(1 - η)s2 - 12η(1 + 2η) + 18η2s (14)
(6)
Equations 7-14 can be used to calculate the two integrals in eq 10. First we carry out the calculations for the integral containing the L-J energy of interaction. Let UL-J(s) be the inverse Laplace transform of r*u*L-J(r*); then by definition,
The perturbation interaction energies can include the effect of various forces, such as dispersion, dipole, quadrupole, and induction forces. Due to the primitive nature of the model employed and the very large dipole moment of amino acids (Greenstein and Winitz, 1961; Khoshkbarchi and Vera, 1996b), only the contributions of dispersion and D-D interactions are considered here. In this work the dispersion forces are represented by the L-J (6-12) model, uL-J(r), written in a dimensionless form, u*L-J(r*), as (Maitland et al., 1981)
u*L-J(r*) )
where D denotes the dipole moment of the component, 0 is the permittivity of the vacuum, and r is the relative dielectric constant of the medium. Assuming the additivity of the energies of interaction which implies
(5)
with
r* ) r/σ
(8)
0 r
APer ) 2πFσ3[ NkT
where σ is the size parameter, F is the number density of the solute, T is the absolute temperature, k is the Boltzmann constant, and superscript hs denotes the hard sphere system. In the primitive framework, the number density corresponds to the number of molecules of solute per unit volume of the system. The contribution of the perturbation term to the residual Helmholtz free energy of the system is represented by an expression based on the Barker-Henderson perturbation theory (Barker and Henderson, 1967) truncated after the first perturbation term. Here, we use the first-order Barker-Henderson perturbation theory for a singlecomponent primitive system as proposed by Tiepel and Gubbins (1973). The equation is written as
APer ) 2πF NkT
uD-D(r) D4σ6 )r*-6 2 kT 3(4π kT)
eq 5 can be written as
with
π η ) Fσ3 6
u*D-D(r*) )
(7)
where is the depth of the potential well according to the L-J model. The D-D interactions are represented by a dimensionless angle-averaged Keesom D-D interactions model as (Maitland et al., 1981)
r*u*L-J(r*) )
∫0∞UL-J(s) e-sr
*
ds
(15)
Multiplying both sides of eq 15 by r*ghs(r*) and integrating between the limits 1 and infinity, after combining with eq 11, we obtain
∫1∞u*L-J(r*) ghs(r*) r*2 dr* ) ∫0∞G(s) UL-J(s) ds
(16)
As eq 16 indicates, once the inverse Laplace transform of r*u*L-J(r*) is known, the corresponding integral in the perturbed term can be easily evaluated. The integral containing the D-D interactions can be similarly evaluated. The inverse Laplace transforms for the L-J and D-D energies of interactions are calculated as
3054 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998
UD-D(s) ) L -1[r*u*D-D(r*)] )
(
) ( )
4 s10 s4 kT 10! 4!
(17)
D4σ6 s10 3(4π0rkT)2 4!
(18)
UL-J(s) ) L -1[r*u*L-R(r*)] )
Combining eqs 1-18 provides us with an expression for the residual Helmholtz free energy of an amino acid in the solution. Using the expression for the residual Helmholtz free energy, the residual chemical potential of the amino acid can be calculated from the following exact thermodynamic relation:
µri )
( )
1 ∂Ar V ∂Fi
(19) V,T
The unsymmetric activity coefficient of a solute, γi, is then related to its residual chemical potential by the following thermodynamic relation (Khoshkbarchi and Vera, 1996b):
ln γi )
µri kT
(20)
One of the advantages of using a primitive model is that, in eq 20, the ideal solution behavior corresponds to the equivalent of the ‘ideal gas’ limit of the hard sphere equation of state. At infinite dilution, the number density of the solute, F, tends to zero, and all attractive energies between the solute molecules are zero. Thus, the residual chemical potential is null, and the activity coefficient of the solute, given by eq 20, is normalized to unity with reference to the state of infinite dilution. It should be mentioned that, in addition to the interaction energies employed in the model proposed here, other forms of the interaction energies can be incorporated in the perturbed term. This provides flexibility for the application of the model to more complex systems, such as those composed of mixtures of amino acids and amino acid-water-electrolyte systems. The formulation just presented is for a single solute; an extension of the model to multicomponent systems is presented in the Appendix. Evaluation of Parameters To correlate the activity coefficients of amino acids and peptides in aqueous solutions using the model proposed in this study, it is necessary to evaluate the size parameters, the depths of potential wells in the L-J energy expression and the dipole moments. The values of the dipole moments of amino acids and peptides were obtained from our previous study (Khoshkbarchi and Vera, 1996b). These values were calculated from a quantum mechanical approach using the Hyperchem Molecular Modeling software. Details of the calculation method are presented elsewhere (Khoshkbarchi and Vera, 1996b). To justify the significance of the parameters obtained from the proposed model, we consider here the GibbsBogoliubov inequality (Hansen and McDonald, 1974):
Ar,hs Ar e + 2πFσ3 NkT NkT
∫1∞u*(r*) ghs(r*) r*2 dr*
(21)
Notably, in its lower bound, this inequality corresponds to the form of the model as given by eqs 1-5.
Figure 1. Activity coefficients of five amino acids at various molalities. Key: (b) experimental data (Fasman, 1976); (-) result of the modeling.
Thus, by equating the left-hand side of the inequality in eq 21 with the experimental value of the residual Helmholtz energy, the adjustable parameters used in the right-hand side should have values that minimize the difference between both sides of the inequality. To achieve this, the size parameter and the depth of the potential well of the amino acids were evaluated minimizing the following objective function:
O.F. ) (γiexp - γical)2
(22)
The number density F was approximated by the molality of the amino acid in the solution. This approximation makes the application of the model simpler because the experimental data for activity coefficients are usually reported in terms of molality. Amino acids, when dissolved in aqueous solutions form ionic species due to the ionization of their carboxyl and amino groups. After dissolving in water, the carboxyl group of the amino acid loses a proton and becomes negatively charged, whereas the amino group gains a proton and becomes positively charged. As a result, an amino acid molecule possesses, at the same time, a negative and a positive charge and therefore carries no net charge. In the presence of a proton donor agent, the amino acid carboxyl group gains a proton and the amino acid molecule becomes positively charged, whereas in the presence of a proton acceptor agent, the amino group of the amino acid loses a proton and the amino acid molecule becomes negatively charged. In the absence of a proton donor or acceptor, up to 99% of the amino acid molecules are in their doubly charged form (Cohn and Edsall, 1965). Therefore, it is reasonable to assume that, in the absence of a proton donor or acceptor, the properties of different ionized forms of the amino acids in aqueous solutions are those of the doubly charged form. Results and Discussion The model developed here was employed to correlate the activity coefficient data of several amino acids and peptides reported in the literature (Fasman, 1976). Figures 1 and 2 show the results of the correlation of the activity coefficients of several amino acids, at their natural pH, as a function of molality. Two different figures are presented due to the different range of molalities required by the data. Figure 3 shows the result of the correlation of the activity coefficients of four simple peptides as a function of their molalities. As
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3055 Table 2. Pooled Standard Deviations (psd) for 95% Confidence Interval of the Correlation of the Activity Coefficients Using the Proposed Model with the Perturbation Parameters Reported psd for perturbation parameters
Figure 2. Activity coefficients of four amino acids at various molalities. Key: (b) experimental data (Fasman, 1976); (-) result of the modeling.
Figure 3. Activity coefficients of four peptides at various molalities. Key: (b) experimental data (Fasman, 1976); (-) result of the modeling. Table 1. Values of σ and E/k for Several Amino Acids and Peptides and Root-Mean-Square Deviation (rmsd) of the Correlation of Experimental Activity Coefficient Data amino acid
/k (K)
σ × 1010 (m)
rmsd × 100
alanine R-amino-n-butyric acid glycine hydroxyproline proline serine threonine valine β-alanine γ-amino-n-butyric acid -amino-n-caproic acid alanylalanine alanylglycine glycylalanine glycylglycine
97.24 98.68 191.11 139.20 12.36 209.87 136.00 42.29 100.25 154.50 101.27 152.10 185.04 180.25 232.15
7.11 7.30 6.61 6.98 4.90 8.01 7.65 7.50 8.60 8.28 8.58 11.92 11.28 11.12 10.30
0.51 0.42 0.84 0.15 7.12 1.19 0.25 0.45 0.31 0.41 0.88 0.40 0.71 0.52 2.81
shown in these figures, the model can accurately represent the activity coefficient experimental data of amino acids and peptides in aqueous solutions. The values of the regressed size parameter, depth of the potential well and the root-mean-square deviation (rmsd) of the correlation of the experimental data for amino acids and peptides studied are presented in Table 1. For proline, an amino acid with large solubility in water, the deviation is probably due to the use of the density of water instead of the density of the solution. A comparison of the calculated rmsd produced by the model developed here with the values obtained using a simplified version proposed previously (Khoshkbarchi
amino acid
/k
σ
glycine alanine serine glycylglycine
(0.004 (0.003 (0.007 (0.008
(0.007 (0.008 (0.013 (0.009
and Vera, 1996b) shows that the simplified model, despite its naive radial distribution function expression, can correlate the experimental data with similar, and in some cases slightly higher, accuracy. This accuracy may be due to the complexity of the present model which in turn renders its numerical optimization difficult. However, the size parameters and depth of potential wells evaluated from the present model, owing to its stronger theoretical basis, are more physically meaningful. To examine the sensitivity of the model to the value of the parameters and the effect of the round off error on the results of the correlations, a sensitivity test for three sample amino acid systems and one peptide system was performed. The results for the activity coefficient obtained from the model were collected for five different runs by perturbing the value of one of the parameters in the interval of the (0.01 while keeping the other parameter constant at its reported value. The pooled standard deviation with 95% confidence interval of the collected results were then calculated. The results of this calculation, reported in Table 2, give an indication of the stability of the present model and show that the round off errors of the parameters do not affect the results. Conclusions A theoretical model based on the perturbation theory of a P-Y hard sphere reference system has been developed to correlate the activity coefficients of several amino acids and peptides in aqueous solutions. The introduction of a more realistic radial distribution function in the perturbation term required more sophisticated mathematics to solve the integrations, but the final expressions are simple to use. Comparison of the model with the Gibbs-Bogoliubov inequality justifies the method used for fitting the parameters. The model accurately correlates the activity coefficients of amino acids and peptides over a wide range of concentration. The results of the model indicate that due to high dipole moment of amino acids and peptides, the D-D interactions play a central role in the chemical potential of amino acids. The comparison of the results obtained with the present model with those given by a simplified version presented previously (Khoshkbarchi and Vera, 1996b) suggests that although the present model is more rigorous, for all engineering applications, the simplified version is adequate. Acknowledgment The authors are grateful to the Thermodynamics Department at Hyprotech for helpful discussions and calculations and to NSERC Canada for financial support. Notation A ) Helmholtz free energy D ) dipole moment
3056 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 det(s) ) dummy function of eq A5, defined by eq A6 g(r) ) radial distribution function i, j, l, m ) component indices k ) Boltzmann constant (1.381 × 10-23 J K-1) N ) number of molecules O.F. ) objective function P ) pressure R ) universal gas constant r ) intermolecular distance s ) Laplace transform operator T ) absolute temperature u ) energy of interaction V ) molar volume
written as
µPer ) 2π i
Greek Letters γ ) activity coefficient ) depth of potential well 0 ) permittivity of vacuum (8.854 × 10-12 C V-1 m-1) r ) relative dielectric constant µ ) chemical potential η ) packing factor F ) number density σ ) size parameter φ1, φ2 ) dummy functions defined by eqs A7 and A8, respectively
kT
) - ln(1 - η3) +
∫0∞rgij(r) e-sr dr )
[( 1
(
)
2π
π2
3η2σi + 3η1σi πP σi + + 6kT 1 - ξ3 (A1)
m
[
(n ) 0, 1, 2, 3)
(A2)
]
η23(3 - η3) η0 3η1η2 6kT + + π (1 - η3) (1 - η )2 (1 - η3)3 3
∑∑n Fmφ1(σm)(σm - σi)(σm - σi)
φ2(σi) )
1 - sσi - e-sσi
(A6)
(A7)
s2
1 - sσi - e-sσi + 0.5s2σi2 s2
(A8)
The interaction energies in the perturbation term can include the effect of dispersion forces, D-D, and dipoleinduced dipole terms, etc. The dispersion forces are represented by a LennardJones (6-12) model (L-J) as (Maitland et al., 1981)
σij12 r12
-
)
σij6 r6
(A9)
where mixing rules are required for σij and ij. The dipole-dipole interactions are represented by an angle-averaged Keesom dipole-dipole interaction model as (Maitland et al., 1981)
(A3)
The contribution of the perturbed term to the residual chemical potential can be represented by a first-order Barker-Henderson perturbation theory (Barker and Henderson, 1967). Following Tiepel and Gubbins (1973), this equation for a multicomponent system can be
+
2(1 - η3)2 m
φ1(σi) )
(A5)
+
2(1 - η3)
4(1 - η3)
(
n
πη3
πη2σm2
uL-J ij (r) ) 4ij
∑ Fkσkn 6k)1
]
[∑ ( ) ( )]
(1 - η3)
Fmφ2(σm) 1 +
2
where
π
+
∑Fmφ1(σm)(σm - σi)(σm - σi)
Fmφ1(σm) σm +
(1 - ξ3)2
ηn )
)
πη3
1+
(1 - η3) det(s) s2 2(1 - η3) σ σ πη 2 i j 1 σij + + s 4(1 - η3)
det(s) ) 1 -
3
(A4)
jl
e-sσij
where
9η22σi2
Phs )
∞
2(1 - η3)s m
The perturbed hard sphere model proposed in this study is generalized to be able to model the activity coefficient of amino acids in multicomponent systems. Due to the complex nature of the equations and to simplify the utilization of them for further calculations, the expressions are directly given for the residual chemical potential. The P-Y residual chemical potential is written as (Reed and Gubbins, 1973) hs
∂Fi
n
∑ ∑FjFl∫σ ujl(r) gjlhs(r) r2 dr] j)1 l)1
{G(s)}ij ) L [rgij(r)] )
π
Appendix
µr,hs i
[
In eq A4, r is the intermolecular distance, µjl(r) is the interaction energy as a function of intermolecular distance, and ghs jl (r) represents the radial distribution function of the hard sphere reference system. The sums run over all components. It should be noted the terms within bracket in eq A4 correspond to the expression for residual Helmholtz free energy. The Laplace transform of rghs jl (r) obtained from the solution of the P-Y equation for hard sphere systems is written as (Blum and Høye, 1977)
Superscripts cal ) calculated D-D ) dipole-dipole interaction term exp ) experimental hs ) hard sphere L-J ) Lennard-Jones Per ) perturbed term r ) residual property Ref ) reference term * ) dimensionless character
n
∂
uijD-D(r) ) -
Di2Dj2 3(4π0r)2kTr6
(A10)
Literature Cited Barker, J. A.; Henderson, D. Perturbation Theory and Equation of State for Fluids. II. A Successful Theory of Fluids. J. Chem. Phys. 1967, 47, 4714.
Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 3057 Blum, L.; Høye, J. S. Mean Spherical Model for Asymmetric Electrolytes. 2. Thermodynamic Properties and the Pair Correlation Function. J. Phys. Chem. 1977, 81, 1311. Bokis, C. P.; Donohue, M. D. Equations of State for Mixtures of Square-Well Molecules: Perturbation and Local Composition Theories for Binary Mixtures of Equal-Size Molecules. Ind. Eng. Res. 1995, 34, 4553. Cohn, E. J.; Edsall, J. T. Proteins, Amino Acids and Peptides as Ions and Dipolar Ions; Hafner Publisher: New York, 1965. Fasman, G. D. Handbook of Biochemistry and Molecular Biology, 3rd ed.; Physical and Chemical Data; CRC: Cleveland, 1976, Vol. 1. Greenstein, J. P.; Winitz, M. Chemistry of Amino Acids, Vol. 1; John Wiley & Sons: New York, 1961. Grundke, E. W.; Henderson, D.; Barker, J. A.; Leonard, P. J. Perturbation Theory and the Excess Properties of the Mixtures of Simple Liquids. Mol. Phys. 1973, 25, 883. Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: New York, 1976. Henderson, D. Perturbation Theory of Mixture of Hard Spheres and Square-Well Molecules. J. Chem. Phys. 1974, 61, 926. Khoshkbarchi, M. K.; Vera, J. H. Measurement and Modelling of Activities of Amino Acids in Aqueous Salt Systems. AIChE J. 1996a, 42, 2354. Khoshkbarchi, M. K.; Vera, J. H. A Simplified Hard-Sphere Model for Activity Coefficients of Amino Acids and Peptides in Aqueous Solutions. Ind. Eng. Chem. Res. 1996b, 35, 4319.
Lebowitz, J. L. Exact Solution of Generalized Percus-Yevick Equation for a Mixture of Hard Spheres. Phys. Rev. 1964, 133, A895. Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular forces; Clarendon: Oxford, 1981. McQuarrie, D. A.; Katz, J. L. High-Temperature Equation of State. J. Chem. Phys. 1966, 44, 2393. Reed, T. M.; Gubbins, K. E. Applied Statistical Mechanics; McGraw-Hill: New York, 1973. Rowlinson, J. S. The Statistical Mechanics of Systems with Steep Intermolecular Potentials. Mol. Phys. 1964, 8, 107. Tiepel, E. W.; Gubbins, K. E. Thermodynamic Properties of Gases Dissolved in Electrolyte Solutions. Ind. Eng. Chem. Fundam. 1973, 12, 18. Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method-I. Nonpolar Gases. J. Chem. Phys. 1954, 22, 1420.
Received for review November 18, 1997 Revised manuscript received February 6, 1998 Accepted February 9, 1998 IE970806G
