17638
J. Phys. Chem. B 2006, 110, 17638-17644
Corresponding-States Laws for Protein Solutions Panagiotis Katsonis,† Simon Brandon,§ and Peter G. Vekilov*,†,‡ Departments of Chemical Engineering and Chemistry, UniVersity of Houston, Houston, Texas 77204, and Department of Chemical Engineering, Technion-IIT, Haifa 32000, Israel ReceiVed: May 2, 2006; In Final Form: July 6, 2006
The solvent around protein molecules in solutions is structured and this structuring introduces a repulsion in the intermolecular interaction potential at intermediate separations. We use Monte Carlo simulations with isotropic, pair-additive systems interacting with such potentials. We test if the liquid-liquid and liquidsolid phase lines in model protein solutions can be predicted from universal curves and a pair of experimentally determined parameters, as done for atomic and colloid materials using several laws of corresponding states. As predictors, we test three properties at the critical point for liquid-liquid separation: temperature, as in the original van der Waals law, the second virial coefficient, and a modified second virial coefficient, all paired with the critical volume fraction. We find that the van der Waals law is best obeyed and appears more general than its original formulation: A single universal curve describes all tested nonconformal isotropic pair-additive systems. Published experimental data for the liquid-liquid equilibrium for several proteins at various conditions follow a single van der Waals curve. For the solid-liquid equilibrium, we find that no single system property serves as its predictor. We go beyond corresponding-states correlations and put forth semiempirical laws, which allow prediction of the critical temperature and volume fraction solely based on the range of attraction of the intermolecular interaction potential.
Introduction cataracts,1
Several human diseases such as eye sickle cell anemia,2 Alzheimer’s,3 the prion diseases,4 and others are associated with the formation of protein condensed phases: liquid or solid. The best studied of the phase transitions in protein solutions is solidification, where protein molecules, instead of being homogeneously distributed in the solution, aggregate into compact structures that can be either ordered (e.g., crystals) or unordered (e.g., amorphous). Another phase transition in protein solutions is liquid-liquid separation, where protein-rich and protein-poor solution phases coexist.5-7 Although the nature of solidification and liquid-liquid separation is far from completely understood, some of their essential features have been systematized, and several general laws of protein phase behavior have been formulated. In certain aspects, these two phase transitions are analogous to the liquid-solid and gas-liquid transitions observed in simple one-component fluids.8 The analogy is similar to the analogy between colloid solutions and gases and liquids, highlighted by Einstein9 and Perrin;10 on the basis of this analogy, protein solutions are often considered in theory as one-component systems consisting of protein particles only. The first model of gas-liquid-phase behavior in onecomponent systems was put forth by van der Waals, who accounted for molecular sizes and interactions in his equation of state.11 Later, he showed that the phase behavior does not depend on the values of the size and attraction parameters if the phase equilibria are described in “reduced” values of temperature T, pressure P, and volume V, i.e., scaled with their * Author to whom correspondence should be addressed. E-mail:
[email protected]. † Department of Chemical Engineering, University of Houston. ‡ Department of Chemistry, University of Houston. § Department of Chemical Engineering, Technion-IIT.
respective values Tcr, Pcr, and Vcr at the critical point for gasliquid coexistence.12 This “law of corresponding states” is an illustration of the power of physics to quantitatively predict phenomena in dramatically different materials; the limitation of the van der Waals formulation is that all systems subject to that law must satisfy a specific equation of state. In further developments, the law of corresponding states acquired a more general foundation. Instead of an equation of state, the modeled systems were characterized by their molecular interactions and kinetic energy: Particles of mass m and hardcore diameter σ interact via a force constant , so that the force f between two molecules at distance r is f(r) ) ψ(r/σ), where ψ is an analytical function. Assuming that temperature is proportional to the average molecular kinetic energy, “corresponding states of motion” were defined for different substances.13,14 This law is more general than the van der Waals law: It applies to systems interacting with any type of isotropic potential. It is restricted to groups of systems characterized by a potential function ψ with varying depth of the attractive minimum, i.e., to so-called conformal systems. In a more recent corresponding-states theory,15,16 the reduced temperature Tr is replaced by the inverse of a new variable X(Tr) ) (1/Tr - Tcr/TB)/(1 - Tcr/TB) where at TB (Boyle temperature) the second virial coefficient B2 becomes zero. According to this definition, 1/X(Tr) is an approximation for the reduced B2/B2cr. This law was found to apply not only to spherical but also to anisotropic and chainlike molecules. This theory is the first in which the second virial coefficient B2 is used instead of a macroscopic thermodynamic variable, and it shows that B2 is a good measure of the coupling between the temperature and the intermolecular potential.15,16 The power of the corresponding-states laws is that a coexistence line for just one of a series of conformal systems could be used for the prediction of all of them provided that their critical
10.1021/jp062698u CCC: $33.50 © 2006 American Chemical Society Published on Web 08/12/2006
Corresponding-States Laws for Protein Solutions
J. Phys. Chem. B, Vol. 110, No. 35, 2006 17639 metric particles as a model for protein molecules, which are often highly anisotropic, makes the conclusions expected from our efforts semiquantitative at best. Still, these conclusions would provide guidance as to how variations in the solution composition, which affect the water structuring and the other factors for protein intermolecular interactions, should affect the protein phase equilibria. Below, we first test if the corresponding-states laws for liquid-liquid coexistence apply to protein-like systems, interacting via potentials with multiple extrema. Then we test if the liquid-solid coexistence in such protein-like systems can be predicted. For these two goals, our main tool is Monte Carlo simulation of the phase diagrams of systems interacting via potentials representative of protein solutions. We also use recent experimental data on the liquid-liquid separation at several compositions of solutions of the protein lysozyme to test their compliance with the corresponding-states law and their correspondence to modeled phase diagrams. Furthermore, we analytically derive and test, using Monte Carlo results, the extension to these laws that allows prediction of the critical properties based on a parameter of the intermolecular interactions. Methods
Figure 1. Potential functions used in the simulations, described by eqs 1 and 2. (a) The “hump” function h(r) ) 0 for all r’s, see eq 2; the value of R is shown below each curve. (b) R ) 50; the values of the parameters γ and rst in the definition of h(R), eq 2, are shown. The thick line represents h(r) ) 0.
properties are known. An extension of the corresponding-states behavior was recently suggested on the basis of Monte Carlo simulations results for several, nonconformal potentials, which allows predictions of the unknown critical temperature and volume fraction.17 The liquid-solid coexistence has proven more difficult to discus in a unified fashion. While for conformal substances the use of temperature and density scaled with the critical values for gas-liquid coexistence results in overlaying liquid-solid coexistence lines, for nonconformal potentials this is not so. Separate laws were proposed for the prediction of liquid-solid equilibrium, the most prominent among them relying on the use of B2 or related parameters.18,19 The objective of the investigations discussed here is to test the applicability of the above laws to systems representative of protein solutions. In protein solutions, especially at ionic strengths higher than 0.1 M, the ranges of intermolecular interactions are significantly shorter than the sizes of the molecules.5,20 For ionic strengths close to physiological (0.100.15 M), the characteristic Debye length κ-1 is between 9 and 6 Å 21 and even shorter at the higher ionic strengths typically used in protein crystallization. Thus, the electrostatic doublelayer interactions are less important for protein molecules in solution, and the shorter-range repulsive maximum, due to structuring of the solvent at the molecular surface, is more important. This structuring is added to the action of a number of interaction forces (van der Waals, double layer).22-26 Such a combination of forces has been shown to lead to nontrivial twobody interaction potentials, involving multiple extrema.22,27,28 These modifications are adequately modeled by the introduction in the interaction potentials of a repulsive maximum at separations from the protein molecular surfaces equal to several water sizes, i.e., longer than the separation for the attractive minimum underlying the formation of condensed phases: crystals, fibers, gels, dense liquids, etc.22,23,28-30 The use of spherically sym-
Model. The interactions between protein molecules in solution were modeled by a generalization of the Lennard-Jones potential for protein interactions that was introduced in ref 31
{{
∞ 1 + h(r) 4 R U(r) ) 2 2 6 R [(r/σ) - 1] [(r/σ)2 - 1]3
}
reσ σ>r
(1)
where the “hump” function
h(r) )
{
r e rst 0 γ[(r/σ)2 - (rst/σ)2]2 r > rst
(2)
introduces a maximum or a secondary minimum, depending on the sign of γ, at intermediate separations. This potential allows independent variation of the range of interactions via the value of R (ref 5) and for the energy barriers imposed by water structured around the protein molecules22,23,28-30 via the function h(r). A negative value of γ in h(r) yields a secondary minimum that can also be attributed to the solvent effects.32 Figure 1 depicts all of the intermolecular potential functions that were used in this work. In Figure 1a only the value of R varies, altering the effective molecular size and the effective range of attraction. In Figure 1b, the values of γ and rst vary, altering the effective range of attraction and the shape of the potential function. The potentials in Figure 1 are not conformal: They do not transform into one another by scaling of their depth . Below, we use dimensionless temperate T ≡ kBT/, where kB is the Boltzmann constant, distance r ≡ r/σ, and energy U ≡ U/. Monte Carlo Simulations. To estimate the liquid-liquid coexistence lines, we use the Gibbs ensemble method introduced in ref 33. Two systems at the same temperature exchange particles and volume, although the total number of particles and the total volume are kept constant. If phase separation is thermodynamically favorable, then at equilibrium each system consists of only one phase (protein-rich or protein-poor): The surface formation energy cost prevents the mixing of the two phases. However, in the vicinity of the critical point, as
17640 J. Phys. Chem. B, Vol. 110, No. 35, 2006
Katsonis et al.
Figure 2. Liquid-liquid coexistence lines for the potentials of Figure 1. (a) h(r) ) 0; the value of R is shown. (b) R ) 50; the values of γ and rst are shown.
temperature increases, this energy cost decreases and vanishes, and the error of the estimated properties becomes larger. In a typical simulation, 300 particles are placed in each system. The simulation length for each phase diagram point varies from 60 000 to 200 000 Monte Carlo steps, where an average step consists of 1 displacement attempt per particle and 3 volumechange and 50 particle-exchange attempts. Systems with shorter ranges of attraction seem to equilibrate slower than systems with longer ranges, so longer simulation times were used. To estimate the liquid-solid boundaries we applied the Gibbs-Duhem integration method introduced in ref 34. Here, two systems, each with a constant number of particles, are at the same temperature and under the same pressure. The pressure is not independent of temperature but is obtained by integration of the Gibbs-Duhem equation. To start the integration, we use an initial point obtained at infinite temperature (T ≈ 1010) where non-interacting hard-sphere behavior is assumed. The initial configurations were 108 particles in a face-centered cubic (FCC) arrangement for the system that represents the solid phase and 128 randomly distributed particles for the system that represents the liquid phase. Typically, the simulation length was 25 000 Monte Carlo steps; an average step consists of 1 displacement attempt per particle and 5 volume-change attempts. Since the accuracy of the simulations in the critical area is low, the critical temperature Tcr and density Fcr were determined from fits of data obtained at temperatures below the critical to the relations33,35
(FL + FG)/2 ) Fcr + A(Tcr - T)
(3)
FL - FG ∝ (Tcr - T)β In eqs 3, A is the critical amplitude, and we used a nonclassical characteristic exponent β ) 0.325.36 Results and Discussion Laws of Corresponding States for Liquid-Liquid Coexistence. Figure 2a shows the liquid-liquid coexistence line as
Figure 3. Liquid-liquid coexistence for all studied potentials. Open symbols, h(r) ) 0; the value of R is shown in plot. Filled symbols, R ) 50; the values of γ and rst are shown. (a) The second virial coefficient B2 scaled with the value of B2 for non-interactive hard spheres of the same effective size as a function of volume fraction φ. (b) Reduced second virial coefficient B2r as a function of reduced volume fraction φr.
it shifts in response to varying the effective particle size and the range of interactions through the value of R in the potential function in eq 1. Figure 2b shows how this line shifts when the shape and the interaction range change due to changes in the hump function h. As expected,8 systems with “broader” attractive ranges have higher critical temperatures and lower critical densities. A measure of the broadness of the attraction or repulsion is the second virial coefficient B2, which depends on temperature T and the intermolecular interaction potential U(r) according to
B2 ) 2π
∫0∞ r2[1 - e-(r)/T] dr
(4)
First, we test if the data in Figure 2 comply with a simple law of corresponding states, formulated in terms of reduced B2 and reduced volume fraction φ. In Figure 3 the temperature axis of the phase diagram has been replaced by the normalized B2n ) B2/B2hs. B2hs characterizes hard-sphere systems with the same effective particle diameter, B2hs ) 2πσeff3/3, where the effective diameter σeff is defined using the collision diameter d as37
σeff ≡
∫0d [1 - e-U(r)/T] dr
(5)
Figure 3a shows that the critical B2 and φ vary significantly: B2 between -1.6B2hs and -1.2B2hs and φcr between 0.18 and 0.36. In Figure 3b, we have scaled the coexistence lines from Figure 3a with respect to these critical values B2cr and φcr. Contrary to observations with a limited number of tested systems,19,38 predictions of the critical temperature and concen-
Corresponding-States Laws for Protein Solutions
J. Phys. Chem. B, Vol. 110, No. 35, 2006 17641
Figure 5. Dependence of the reduced second virial coefficient Br ) B2/B2cr on reduced temperature Tr for all studied potential functions. For symbol explanation, see Figures 2 and 3.
Figure 4. Laws of corresponding states. (a) van der Waals formulation. The line is just a guide for the eye. (b) Hall and Hacker formulation; see text. For symbol explanation, see Figures 2 and 3.
tration for liquid-liquid coexistence or the liquid-liquid coexistence line, based on B2, are of limited accuracy.39 To test how efficiently the van der Waals law of corresponding states applies to the studied cases, the liquid-liquid coexistence data were plotted in Tr-φr coordinates in Figure 4a: A universal curve emerges. Comparing this curve with published liquid-liquid coexistence lines for the simplest attractive potential model, the square-well potential
Usw(r) )
{
∞ r
