Theoretical study of phase behaviour of DLVO model for lysozyme and $\gamma$-crystalline aqueous electrolyte solutions

Mean spherical approximation (MSA), second-order Barker-Henderson (BH) perturbation theory and thermodynamic perturbation theory (TPT) for associating fluids in combination with BH perturbation theory are applied to the study of the structural properties and phase behaviour of the Derjaguin-Landau-Verwey-Overbeek (DLVO) model of lysozyme and $\gamma$-cristalline aqueous electrolyte solutions. Predictions of the MSA for the structure factors are in good agreement with the corresponding computer simulation predictions. The agreement between theoretical results for the liquid-gas phase diagram and the corresponding results of the experiment and computer simulation is less satisfactory, with predictions of the combined BH-TPT approach being the most accurate.


Introduction
Among globular proteins, the equilibrium properties (structure factors, thermodynamic properties and phase behavior) of aqueus solutions of lysozyme and γ-crystalline are the ones, perhaps, most thoroughly studied. Numerous recent studies have been reported for these systems from both theoretical and experimental perspective (see [1][2][3][4][5][6] and references therein). Much effort, in particular, has focused on the investigation of the phase behavior of the lysozyme solution. Experimentally, the liquid-gas phase coexistence for this system was described by Ishimoto and Tanaka [7]. According to the later studies of Broide et al. [8], this phase coexistence is unstable with respect to crystallization. Worth mentioning is an important contribution due to George and Wilson [9], who discovered the existence of the 'crystallization slot' for the values of the second osmotic virial coefficient in the vicinity of the liquid-gas critical point, where one might expect crystallization of the proteins.
Since globular proteins can be viewed as colloidal macroions, most of the theoretical studies of protein phase equilibrium are based on the concepts borrowed from the physics of colloids. As far as the interaction between two protein macromolecules is very complicated and, to a large degree, the known theoretical studies are based on the coarse-grained potential models [10]. The simplest version of the model, the so-called one-component model, represents the effect of the solvent and electrolyte produced by the continuum approximation. Usually, the corresponding effective interaction is represented by the hard-sphere (or soft-sphere) interaction combined with long-range repulsive screened Coulomb interaction and short-range attractive van der Waals interaction. This is the model utilized in the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory of the colloidal stability [11]. Recent computer simulation studies [4,5,12,13] demonstrate that DLVO model, being not very accurate in reproducing the short-range behavior of the experimental structure factors at higher pH values, appears to be quite successful in describing the phase behavior of the lysozyme solution. In particular, DLVO model was capable of reproducing the flat portion of the experimental phase diagram in the vicinity of the critical point. Further progress in the coarse-grained modelling of protein solutions is related to the introduction of the orientationally dependent short-range attraction between colloidal particles [14][15][16]. The models of this type are aimed at a more detailed description of the protein molecules taking into account the existence of charged groups on their surfaces. Phase behavior and structural properties of these models have been studied in references [3][4][5][17][18][19][20][21][22][23][24]. Further steps in the detalization of the interaction in protein solutions are connected with the substitution of the one-component model with the one that contains multiple components and takes into account simple salt ions and solvent molecules. The presence of a simple electrolyte can be addressed within the framework of a highly asymmetric electrolyte solution model, which has been extensively studied using integral equation methods (see [25] and references therein) and by explicitly taking into account the effects of association [26][27][28][29][30][31][32]. Finally, the effects due to solvent molecules can be considered by extending the scheme developed for simple electrolytes and polyelectrolytes [33][34][35][36][37][38][39][40] Despite a good performance of the DLVO model in predicting the phase behavior of the electrolyte solution of lysozyme [4,5], the authors have not proposed any theoretical description for its equilibrium properties. The goal of the present study is to fill this gap and develop simple theories capable of describing the structural properties and phase behavior of the model, at least at a qualitative level. With this goal in mind, we apply several simple liquid state theories, including the mean spherical approximation (MSA), second order version of the perturbation theory (TPT), for associating fluids in combination with BH perturbation theory, and critically assess their performance. The paper is organized as follows. In the section 2 we formulate a potential model while in section 3 we discuss the details of the MSA, BH perturbation theory and TPT. Our results are presented in section 4 and our conclusions are given in section 5.

The model
DLVO model utilized in [4,5,11] treats lysozyme solution as an effective one-component fluid of spherical particles with the number density ρ interacting via a pairwise additive potential which is only a function of the interparticle distance r . Pair interaction between particles consist of (i) a short-range attractive part of van der Waals term, where σ DLVO is an effective hard-core diameter and A H is the Hamaker constant, and (ii) a Debye-Hückel where Z e is a net charge of the lysozyme macromolecule in electronic units, 0 is permittivity of vacuum, [5]. The Debye-Hückel screening length κ is defined by the expression where I is the ionic strength of the solution, which takes into account the presence of ions due to the buffer and the added salts. In the present study, we neglect a weak dependence of DH potential on the temperature and assume for κ and r the values calculated at ambient conditions.
The total DLVO potential is written as follows: where the cut-off value δ is introduced to avoid a singularity of the van der Waals contribution at r = σ DLVO and corresponds to the thickness of the Stern layer.

Theory
The properties of the model at hand are studied using MSA, BH perturbation theory and TPT for associating fluids. MSA is a simple analytical approach, which is known for being capable of providing both structural and thermodynamic properties of a large number of model systems with sufficient accuracy. Relatively simple and possibly accurate description of thermodynamic properties of the model can be also achieved within the framework of the BH perturbation theory. Finally, to account for strong and short-range attraction between particles, which characterizes DLVO potential model (4), we apply an appropriate combination of the BH perturbation theory and TPT for associating fluids. The accuracy of each of these approaches is evaluated via comparison of the theoretical predictions with the corresponding computer simulation and experimental predictions.

MSA
MSA consists of the Ornstein-Zernike (OZ) equation whereĉ(k) andĥ(k) are Fourier transforms of the direct and total correlation functions c(r ) and h(r ), respectively, and MSA closure relation For the sake of analytical description, we approximate the vdW part of the potential V vdW (r ) using one-Yukawa potential, i.e., r .
As a result, the pair potential outside the hard core is represented by the two-Yukawa potential and we have used an analytical solution of the MSA due to Blum and Hoye [41]. The structural and thermodynamic properties have been calculated utilizing this solution and using a closed form of analytical expressions presented in [42].

Barker-Henderson perturbation theory
Here, we utilize the second-order BH perturbation theory. Within the framework of BH perturbation theory [43], Helmholtz free energy F of the system per particle is given by the following expression where k B is the Boltzmann constant, F ideal is the ideal gas free energy Λ is the thermal de Broglie wavelength, F HS is the hard-sphere Helmholtz free energy and F 1 and F 2 are the first-and second-order contributions to Helmholtz free energy of the system. Here, for F HS we have used the Carnahan-Starling extension

13604-3
where η = π 6 ρσ 3 0 , σ 0 = σ DLVO + δ and for F 1 and F 2 , we have where κ HS is the compressibility of the hard-sphere fluid. Here, we use the corresponding Carnahan-Starling expression [44]. Using the above expression (8) for Helmholtz free energy of the system, all the rest thermodynamical properties (pressure and chemical potential) can be derived using the standard thermodynamical relations. As for any perturbation theory, the BH model works best if the attractive potential is not large in magnitude. For potentials with strong attractions, other approaches are more suitable.

Thermodynamic perturbation theory for associating fluid
To account for the strong attraction seen in the DLVO potential of lysozyme, we combine the BH perturbation theory and thermodynamic perturbation theory (TPT) for associating fluids with spherically symmetric interaction [45]. Following the earlier studies [26,29,46], the total pair potential of the system (4) is represented as a sum of the reference and associating pieces, i.e., Here, V DLVO (σ 0 ) < V (as) 0 < 0, and we assume that V (as) 0 is temperature dependent. A particular choice for the potential splitting parameter V (as) 0 is discussed below. According to the splitting of the total DLVO potential (13), the system Helmholtz free energy within the framework of the TPT [45] is as follows: where for the free energy of the reference system F ref we have used the second-order expression (8) with DLVO pair potential V DLVO (r ) substituted by the reference potential V (ref) DLVO (r ), and thus for the associative part F (as) we have F as where χ 0 and χ 1 satisfy the following set of equations where

13604-4
Phase behaviour of DLVO model for lysozyme and γ-crystalline aqueous electrolyte solutions m is a maximum number of bonds per particle allowed and g (ref) (σ 0 ) is the contact value of the radial distribution function of the reference system. The latter quantity was obtained using the parametrization of the contact value of the radial distribution function for the hard-sphere square-well fluid [47]. The knowledge of the free energy of a system (15) enables one to calculate thermodynamical properties of interest using standard thermodynamical relations. The version of the TPT approach discussed above has two input parameters, i.e., the maximum number of bonds m and the potential splitting parameter V (as) 0 . In general, the number of nearest neighbours for the model at hand could be up to 12, thus one can assume m = 12. In this case, the value of V (as) 0 , which ensures the saturation of the associative potential V DLVO (r ) for m = 12, can be chosen to be relatively large in comparison with the contact value of the DLVO potential V DLVO (σ 0 ). However, according to the earlier studies [48], 'single bond' approximation utilized here is accurate only for relatively small values of m. In addition, the probability of bonding of 12 particles simultaneously at the densities where the liquid-gas separation occurs is small; thus, the optimal choice for m and V (as) 0 requires a certain compromise. In the present study we assume that m = 3 and the reduced value of the associative potential at the contact V * (as) In addition, to provide an accurate description of the reference system at the critical point, we also assume that V (as) 0 = k B T cr , where T cr is the critical temperature. Combining the latter two assumptions, we For T = T cr we have V (as) 0 = −k B T cr , i.e., for the critical temperature, the minimum value of the reference system potential V (ref) DLVO (r ) is equal to −k B T cr . For this value of the potential minimum, the BH approach is expected to be sufficiently accurate. Although our choice for m and V (as) 0 is rather empirical, the theory proposed is self-contained because there is no need in the input from outside. In particular, critical temperature T cr (and critical density ρ cr ), which enter the expression for V (as) 0 (20), is obtained as usual from the solution of the set of two equations, which requires the first and the second derivatives of the pressure with respect to the density to be equal to zero, i.e.,

Results and discussion
In this section we present our numerical results for the structural properties and phase behavior of lysozyme and γ-crystalline aqueous electrolyte solutions. In both cases, we use the DLVO model with σ DLVO = 37.08 Å for lysozyme and σ DLVO = 37.8 Å for γ-crystalline. For the Stern-layer thickness and Hamaker constant, the following values have been used [5]: δ = 1.8 Å and A H = 18.8 kJ/mol. The structural properties of the lysozyme solution [static structure factor S(k) and radial distribution function g (r )] were calculated using MSA (6). Parameters for the Yukawa potential (7) representing V vdW (r ) were chosen using 'best eye fit' method supplemented by the equality of both potentials at the contact distance, i.e., V Y (σ 0 ) = V vdW (σ 0 ). In figure 1 we compare V Y (r ) and V vdW (r ) with the following choice of the Yukawa potential parameters: A Y = 375.19 kJ/mol·Å and κ Y = 0.45 Å −1 . In figure 2 we compare the corresponding theoretical and computer simulation [5] results for the structure factor S(k) at the values of pH used in the corresponding experimental study, i.e., pH = 2.8, 4.2, 5.07. These values of pH correspond to the following three values of the ionic strength (see table 1 of reference [5]), i.e., I = 0.035 M, 0.081 M, 0.102 M, respectively. In addition, we also show the corresponding experimental results [5] for S(k). Very good agreement is observed between theoretical and computer simulation results for the structure factor at all values of pH and I studied. Although DLVO model is not very accurate in describing the short-range behavior of the structure factor at higher values of pH (figure 2), its predictions for the phase behavior are in reasonable agreement with experimental predictions (see figure 3).  . Experimental values for the critical concentration and temperature for lysozyme [49] and γ-crystalline solutions [50] are (230 ± 10 mg/ml, 273 K) and (289 ± 20 mg/ml, 311 K), respectively. correspond to the protein overall charge Z = 1) and ionic strength I = 0.24 M (red color). These results were obtained using second-order BH theory (8) and our version of the TPT (15). Unfortunately, for the employed parameters of the DLVO model, the MSA does not have a convergent solution in the range of temperatures and densities, where one would expect the location of binodals. In addition, in the same figure we show computer simulation results and the results of the experiment [49,50]. Here, computer simulation predictions are much more accurate in comparison with the predictions for the structure factors demonstrated in the previous figure. Our results obtained using the BH perturbation theory, are much less accurate, giving the values for both critical temperature and critical density that are too large. In addition, BH phase diagram is too narrow in comparison with the computer simulation and experimental phase diagrams. Predictions of our combined BH and thermodynamic perturbation theory are in reasonable agreement with the computer simulation data. Good agreement is observed for the critical temperature. Slightly less accurate are the predictions for the critical concentration. The overall shape of the phase diagram is still too narrow, although it is now much closer to the shape of computer simulation phase diagram.

Conclusions
In this article we studied the structural properties and phase behaviour of the DLVO model of lysozyme and γ-crystalline aqueous electrolyte solutions using MSA, second-order BH perturbation theory and a combined approach based on the BH theory and TPT for associating fluids. Theoretical results are compared with computer simulations and experimental results. Predictions of the MSA for the structure factor of lysozyme solution are in good agreement with the corresponding computer simulation predictions. However, MSA does not have a convergent solution in the range of the temperatures and densities, where one would expect the location of the corresponding MSA binodals. We conclude that MSA is inappropriate for the phase behavior of the DLVO-type models of lyzosyme and, perhaps, other globular proteins. Among the theories used to describe the phase behaviour of the lysozyme and γ-crystalline solutions, only a combined BH-TPT approach provides a reasonable qualitative agreement with computer simulations and experimental description. We expect that a further improvement of the theory can be achieved using association concepts in combination with a more detailed description of the pro-

13604-7
tein molecules taking into account the existence of the charged groups on their surfaces. Corresponding studies are underway and results will be reported in due course [51].