Modeling Opposite Effects of an Additive on Liquid-Liquid Phase Separation and Crystal Solubility of Protein Solutions

1. Introduction

Understanding and controlling condensation of globular proteins in aqueous media is fundamental for comprehending cell compartmentalization,[1,2,3] protein-aggregation diseases,[4,5] developing stable pharmaceutical formulations,[6,7] preparing protein-based materials,[8,9] and generating protein crystals[10,11,12] with applications in structural biology[13] and separation science[14,15].

An interesting phenomenon of protein aqueous mixtures is the reversible formation of metastable protein-rich liquid microdroplets through liquid-liquid phase separation (LLPS), typically induced by lowering temperature.[16,17,18] Since LLPS is a metastable phase transition, the protein-rich phase can act as an intermediate for the formation of other protein condensed phases such as crystals[17,19,20] and aggregates[4,21,22]. This phenomenon was initially considered for just few protein cases: eye-lens crystallins[16,23,24] and lysozyme[17,25]. However, it has recently attracted considerably more attention as it is believed to drive the formation of membraneless organelles in the cytosol and can promote formation of pathological protein aggregates.[1,2,3,4,5] It has also been investigated in the contest of protein pharmaceutical formulations, where LLPS is described to negatively impact the stability and efficacy of monoclonal antibodies.[8] Finally, LLPS is also known to enhance the rate of protein crystallization,[19,26,27,28] which is beneficial not only for the characterization of protein 3D structure[29,30] but also for protein purification in downstream processing[14,15].

From a molecular point of view, all protein condensation processes are driven by solvent-mediated protein-protein attractive interactions.[31] In principle, it is extremely difficult to describe protein-protein interactions as they ultimately depend on the surface composition and spatial orientation of amino-acid groups. These are responsible for the formation of a highly heterogenous distribution of charged chemical groups, hydrophilic moieties and hydrophobic patches. Moreover, additives, which are invariably present in protein solutions, can significantly alter protein-protein interactions.[32] These cosolvents modify these interactions through mechanisms such as electrostatic screening,[33,34] preferential hydration or binding,[35,36] and molecular crowding[18,37,38,39,40,41]. Nevertheless, some key features associated with the collective behavior of globular proteins can be theoretically described by employing a simple one-component colloid model using hard spheres as a reference system.[42,43,44] For example, LLPS metastability with respect to protein crystallization is successfully explained by assuming that the range of protein-protein interactions is relatively short compared to particle diameter.[42]

From a thermodynamic point of view, the phase behavior of a protein aqueous system is typically described by employing a temperature-composition phase diagram. As schematized in Figure 1A, a typical phase diagram shows a protein-crystal solubility curve with the LLPS boundary positioned in the crystal supersaturated domain, at relatively high protein concentrations. Note that the LLPS boundary can be experimentally characterized because protein crystallization kinetics is usually slow.[6,16,42,45]

Additives such as salts, buffer components and polyethylene glycols are typically present in protein aqueous solutions. They are used to stabilize proteins against unfolding, favor protein crystallization and mimic physiological conditions. The thermodynamic effect of additive concentration has often been described using the theory of preferential interaction.[46] Accordingly, an additive that is preferentially excluded from the vicinity of protein surface (preferential hydration) reduces protein solubility (e.g., NaCl for lysozyme, sulfates). In contrast, an additive that accumulates near protein surface (preferential binding) increases protein solubility (e.g., thiocyanates, urea).[47,48] A one-component colloid model may be still used to describe the phase behavior of these multicomponent protein solutions provided that the additive has a low molecular weight and is therefore considered as integral part of the solvent. The effect of additive concentration is then implicitly taken into account by considering a corresponding change in the solvent-mediated protein-protein interactions.[44] For example, attraction energy should increase with additive concentration in the presence of preferential hydration.[49] As illustrated in Figure 1B, this salting-out effect leads to a shift of both crystal solubility and LLPS boundary towards higher temperatures, with crystal solubility shifting toward lower protein concentrations. This behavior has been experimentally demonstrated for lysozyme in the presence of NaCl and other additives.[17,50,51] Interestingly, shifts in the LLPS boundary appear to be stronger than those of crystal solubility. This results in an increase in the metastability gap between LLPS boundary and solubility curve.[51]

In our previous studies, we have shown that addition of 4-(2-hydroxyethyl)-1-piperazineethanesulfonate (HEPES) to lysozyme-NaCl solutions at constant ionic strength and pH produces even a more complex effect on the phase diagram.[52,53] As qualitatively schematized in Figure 1C, HEPES moves the solubility curve toward higher temperatures (as in Figure 1B) while moving the LLPS boundary in the opposite direction. This leads to a significant increase in the metastability gap between the two phase boundaries. Remarkably, addition of HEPES was also found to significantly increase yield of protein crystallization, especially in the presence of LLPS.[28] Clearly, understanding the multifaceted effect of this type of additives the phase behavior of protein solutions, especially in relation to the metastability gap, is important for controlling protein crystallization and other aggregation processes.

To explain the effect of HEPES on lysozyme phase diagram, partitioning experiments were also carried out. It was found that HEPES accumulates in the protein-rich liquid phase showing that HEPES preferentially binds to lysozyme.[52] Furthermore, dynamic-light-scattering experiments showed that HEPES weakens protein-protein attraction energy.[52] While these results are consistent with HEPES suppressing lysozyme LLPS, they cannot explain HEPES salting-out effect on lysozyme solubility.

The lysozyme 3D structure obtained from crystals grown in HEPES buffer also shows that his organic molecules act as a ligand. Indeed, it is found in the lysozyme catalytic site, with the hydroxyethyl and sulfonate groups pointing inward and outward, respectively.[54] The sulfonate group is close to a cationic arginine group of a neighboring protein. It is therefore reasonable to assume that the electrostatic attraction between these two ionic groups may strengthen one of the crystal contacts.[53] In other words, the effect of HEPES on lysozyme solubility may be qualitatively explained by considering that this organic molecule acts as a physical cross-linker enhancing attractive interactions between neighboring protein units in the crystal lattice. This thermodynamically stabilizes the crystalline phase. Our hypothesis aligns with the general belief that multi-functional organic molecules are beneficial in protein crystallography.[55]

In this paper, we examine the ability of a colloid model to quantitatively describe the multifaceted effects of HEPES on lysozyme phase behavior assuming that it does not alter protein conformational state. We will specifically model the effect of this additive by assuming that it causes an increase in protein-protein attractive interactions accompanied with an increase in their anisotropic character. It is also assumed that the additive does not alter protein conformational state.

2. Methods

A protein aqueous solution is described as a one-component compressible fluid of colloidal particles provided that this aqueous medium is regarded as incompressible and low-molecular weight additives such as salts and buffer components are implicitly described as a part of the background solvent. These additives are modulators of solvent-mediated protein-protein interactions. Correspondingly, the protein osmotic pressure, $\mathsf{\Pi }$, becomes a “gas pressure” and LLPS is equivalent to a “gas-liquid” phase transition. Furthermore, protein chemical potential, $\mu$, must describe particle insertion accompanied with isochoric removal of solvent.[31] It is then appropriate to start from the Helmholtz free energy, $A$, of the colloid fluid. It is also practically convenient to introduce the corresponding reduced unitless quantity, $a\equiv \beta (V_{\mathrm{P}}/V)A$, where $\beta \equiv 1/k_{\mathrm{B}}T$, $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ the absolute temperature, $V_{\mathrm{P}}$ the protein volume, and $V$ is the volume of the fluid phase. We shall set $k_{\mathrm{B}}=1$ so that ${\beta }^{\hspace{0.17em}-1}$ is the same as temperature, and energy parameters adopt temperature units. The composition of the fluid phase is described by the volume fraction, $\varphi =N{V}_{\mathrm{P}}/V$, where $N$ is the number of colloidal particles. An experimental volume fraction of lysozyme solution is readily calculated by multiplying lysozyme mass concentration by the known[56] specific volume of this protein (0.713 cm³·g^-1).

The protein phase diagram[42] shows protein volume fraction ($\varphi$) and temperature (${\beta }^{-1}$) on the x- and y-axis, respectively. The LLPS boundary is described by the protein volume fractions, ${\varphi }^{(\mathrm{I})}(\beta )$ and ${\varphi }^{(\mathrm{II})}(\beta )$, of the two coexisting of protein-diluted (“gas”, I) and -concentrated (“liquid”, II) fluid phases. These two volume fractions become equal at the critical point, which is described by the critical coordinates, $({\varphi }_{\mathrm{c}},{\beta }_c)$.To determine ${\varphi }^{(\mathrm{I})}(\beta )$ and ${\varphi }^{(\mathrm{II})}(\beta )$, the expressions of chemical potential and osmotic pressure are extracted from the Helmholtz free energy using

$$\widehat{\mu }={\left({\displaystyle \frac{\partial a}{\partial \varphi }}\right)}_{\beta }$$

(1a)

$$\widehat{\pi }=\varphi \hspace{0.17em}\widehat{\mu }-a$$

(1b)

where $\widehat{\mu }\equiv \beta \mu$ and $\widehat{\pi }\equiv \beta \hspace{0.17em}\mathsf{\Pi }{V}_{\mathrm{P}}$. At any sufficiently high value of $\beta$, the protein volume fractions, ${\varphi }^{(\mathrm{I})}$ and ${\varphi }^{(\mathrm{II})}$, of the two coexisting phases are determined by numerically solving the chemical-equilibrium conditions:

$$\widehat{\mu }({\varphi }^{(\mathrm{I})},\beta )\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\widehat{\mu }({\varphi }^{(\mathrm{II})},\beta )\hspace{0.17em}$$

(2a)

$$\widehat{\pi }({\varphi }^{(\mathrm{I})},\beta )\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\widehat{\pi }({\varphi }^{(\mathrm{II})},\beta )\hspace{0.17em}$$

(2b)

Critical coordinates, $({\varphi }_{\mathrm{c}},{\beta }_{\mathrm{c}})$, are extracted using the conditions: ${(\partial \widehat{\mu }/\partial \varphi )}_{{\beta }_{\mathrm{c}}}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}{({\partial }^2\widehat{\mu }/\partial {\varphi }^2)}_{{\beta }_{\mathrm{c}}}=0$.

The hard-sphere model is employed as a reference model. The reduced Helmholtz free energy, $a$, is then given by

$$a=a_{\mathrm{HS}}+a_{\mathrm{R}}$$

(3)

where $a_{\mathrm{HS}}$ is the hard-sphere contribution and $a_{\mathrm{R}}$ is a residual term describing deviation from reference model. The hard-sphere term, which describes both particle translational motion (ideal contribution) and hard-sphere excluded-volume interactions (steric repulsion), is given by

$$a_{\mathrm{HS}}(\varphi )=\varphi \hspace{0.17em}{\widehat{\mu }}_0+\varphi (\ln \varphi -1)+{\displaystyle \frac{4-3\varphi }{{(1-\varphi )}^2}}{\varphi }^2$$

(4)

where ${\widehat{\mu }}_0$ is the uninfluential standard chemical potential, the second term represents the ideal-gas contribution and the third term describes steric repulsion according to Carnahan-Starling equation of state.[57] It follows from Eqs. 1a,b that the hard-sphere chemical potential and pressure are given by ${\widehat{\mu }}_{\mathrm{HS}}=\hspace{0.17em}{\widehat{\mu }}_0+\ln \varphi +(8-9\varphi +3{\varphi }^2)\varphi /{(1-\varphi )}^3$ and ${\widehat{\pi }}_{\mathrm{HS}}=\varphi (1+\varphi +{\varphi }^2-{\varphi }^3)/{(1-\varphi )}^3$, respectively. The residual term in Eq. 3 characterizes protein-protein attraction energy and vanishes in the limit of high temperatures ($\beta \to 0$).[58] This residual term will be further discussed in the next section.

The crystalline ordered phase, which has a very high protein volume fraction (typically more than 50%), is assumed to be incompressible.[42,59] In this case, the protein chemical potential in the solid crystalline phase, ${\mu }_S(\beta )$, coincides with the Helmholtz free energy of one particle. The crystal solubility boundary is described by the protein volume fraction of the coexisting fluid phase, ${\varphi }_S(\beta )$. At any given $\beta$, this is determined by numerically solving the chemical-equilibrium condition:

$$\widehat{\mu }({\varphi }_S,\beta )\hspace{0.17em}={\widehat{\mu }}_S(\beta )$$

(5)

where ${\widehat{\mu }}_S\equiv \beta {\mu }_S$. A simple expression for ${\widehat{\mu }}_S$ was obtained from cell theory as discussed in the next section.

3. Results and Discussion

To examine the effect of HEPES on the phase diagram of lysozyme aqueous solutions, we have previously considered two systems. The first system contains HEPES (0.10 M) and NaCl (0.15 M) at pH 7.4, while the second system is a reference system (denoted as REF) and mostly contains NaCl (0.183 M; Tris buffer, 0.020 M) at the same pH and ionic strength (0.20 M).[52,53] In this way, replacement of NaCl with HEPES is carried out in conditions in which long-range electrostatic repulsive interactions between proteins remain the same. Note the presence of Tris buffer is REF, which was employed to improve buffer capacity. Although its contribution to the total ionic strength of 0.20 was taken into account, its concentration is fairly low and Tris-specific effects are not expected to be significant.

In the following subsections, we shall first examine the effect of HEPES on lysozyme crystal solubility and then the corresponding effect on the LLPS boundary.

3.1. Crystal Solubility and Thermodynamics of Protein Crystal

To theoretically describe crystal solubility, expressions for the protein chemical potential in the crystal and fluid phases in Eq. 5 are needed. According to cell theory,[42,59] ${\widehat{\mu }}_S$, may be written as:

$${\widehat{\mu }}_S={\widehat{\mu }}_0-n_S{\displaystyle \frac{\beta {\epsilon }_S}{2}}-\ln {\mathsf{\Omega }}_S$$

(6)

where ${\widehat{\mu }}_0$ is the same as in Eq. 4. The second term represents crystal internal energy, with $n_S$ loosely representing the number of interactions “sites” or “contacts” between a central protein and its neighboring surrounding proteins. The factor “2” considers that a protein-protein interaction involves two proteins; i.e., the range of interactions is sufficiently short that one interaction cannot extend to more than two proteins. In Eq. 6, ${\epsilon }_S$ represents the average attraction energy of a contact. In principle, the number of contacts could be defined as involving interactions between single amino acids (atomic contacts) or as relatively big regions on the protein surface encompassing multiple functional groups on the protein surface (region contacts). Thus, the value of ${\epsilon }_S$ depends on the choice of $n_S$ and only their product, $n_S{\epsilon }_S$, can be unambiguously defined. The last term in Eq. 6 is an entropic contribution related to both the residual translational and orientational motion of a protein in the crystal lattice, with ${\mathsf{\Omega }}_S$ representing the phase volume[59] (as a multiple of particle volume, $V_{\mathrm{P}}$) accessible to both particle center of mass and particle rotation. It will be further discussed at the end of this subsection.

The highest experimental[53] ${\varphi }_S$ are 0.014 and 0.027 for the HEPES and REF system, respectively. We assume that these values are sufficiently low that the protein chemical potential of the fluid phase is approximated by its ideal-dilute approximation, $\widehat{\mu }\approx {\widehat{\mu }}_0+\ln \varphi$, neglecting contribution of protein-protein interactions in the fluid phase. We will later validate this assumption by verifying that inclusion of steric repulsion and protein-protein attraction energy terms do not significantly alter solubility boundary. Insertion of Eq. 6 into Eq. 5 yields the following Van’t Hoff equation:

$$\ln {\varphi }_S=-n_S{\displaystyle \frac{\beta \hspace{0.17em}{\epsilon }_S}{2}}-\ln {\mathsf{\Omega }}_S$$

(7)

Van’t Hoff plots for lysozyme solubility in the HEPES and REF systems are shown in Figure 2. Here, we can see that solubility values in the HEPES system are appreciably lower than those in the REF system. Application of the method of least squares yields: $n_S{\epsilon }_S=$(10.4±0.5)10³ K (86.6 kJ·mol^-1) and $\ln {\mathsf{\Omega }}_S=$-(12.9±1.0) for the HEPES system and $n_S{\epsilon }_S=$(10.1±1.0)10³ K (84.0 kJ·mol^-1) and $\ln {\mathsf{\Omega }}_S=$-(12.9±1.8) for the REF system. These values are essentially the same within the experimental errors due to correlation between slope and intercept parameters. To further examine these solubility data, we assume that addition of HEPES has a negligible effect on crystal entropy. This assumption is motivated by the hypothesis that HEPES energetically stabilizes region contacts through physical crosslinking while not appreciably altering crystal lattice, consistent with available crystal structures.[54,60,61] Accordingly, we define $n_S$ as a region contact and assume it is the same for both systems. We then set $\ln {\mathsf{\Omega }}_S=$-12.9 so that we can attribute solubility variations entirely to ${\epsilon }_S$. We can now determine appreciably different values of $n_S{\epsilon }_S$. These are (10.39±0.02)×10³ K and (10.12±0.03)×10³ K for the HEPES and REF systems, respectively. The corresponding linear fits are also shown in Figure 2. These share the same intercept at $\beta$=0.

To determine ${\epsilon }_S$, we need to identify a reasonable value of $n_S$. It is known[61,62] that there are eight lysozyme molecules interacting with every lysozyme molecule in the tetragonal crystal structure. Correspondingly, there are eight region contacts on each protein involved in the formation of multiple interatomic bonds with neighboring proteins. However, the interactions occurring with the two proteins immediately above and below the reference protein (along the crystallographic c axis) are weak and can be neglected.[62] The proteins involved in the interactions with a reference protein (M) are represented on the (0,0,1) crystallographic plane in Figure 3. Here, we can see that there are two neighboring proteins, F (at $-z+1$) and G (at $-z+1/2$), related to M by twofold screw symmetry. There are also two B (at $z+3/4$ and $z-1/4$) and two C (at $z+1/4$ and $z-3/4$) proteins with different positions along the c-axis interacting with M. These are related by fourfold screw symmetry to M.[61,62] The catalytic site (cleft) of protein M, which can be occupied by HEPES, is near the protein B at $z-1/4$. It is therefore reasonable to attribute the increase in attraction energy caused by HEPES to this specific contact region.

In summary, there are six main region contacts on the protein that can be identified as binding sites. We therefore set $n_S=$6 and determine that the average attraction energy of an interacting region is ${\epsilon }_S=$1732 K and 1687 K for the HEPES and REF system, respectively. In summary, our data analysis shows that HEPES causes an increase of 2.7% in the value of ${\epsilon }_S$ inside the crystal.

We conclude this subsection by examining the extracted $\ln {\mathsf{\Omega }}_S$ describing crystallization entropy. Although its accurate interpretation should factor in protein shape and changes in solvent entropy, it is important to examine whether the value of phase volume, ${\mathsf{\Omega }}_S=2.5\times {10}^{-6}$, yields physically acceptable geometric parameters. We specifically assume that ${\mathsf{\Omega }}_S$ is the product of a translational factor, ${\mathsf{\Omega }}_r$, and an orientational factor, ${\mathsf{\Omega }}_{\theta }$. The translational factor is given by ${\mathsf{\Omega }}_r=V_r/V_{\mathrm{P}}$, where $V_r$ is the center-of-mass excursion volume. If proteins are assumed to be spherical particles, we can use unit cell volume of tetragonal lysozyme (238 nm³), number of proteins inside unit cell (8)[61] and protein volume (16.9 nm³) to determine that protein volume fraction in the crystal is 0.57. For spherical particles, translational motion is lost when volume fraction reaches the classical close-packing value of 0.7405. Since an increase of 1.091 in particle diameter is needed to increase experimental volume fraction to close-packing value, the diameter of the excursion sphere is estimate to be $2\times 0.091\hspace{0.17em}\sigma$, where $\sigma$ is the diameter of the protein. This implies that ${\mathsf{\Omega }}_r\approx 6.0\times {10}^{-3}$ from which we calculate: ${\mathsf{\Omega }}_{\theta }\approx 4.1\times {10}^{-4}$. To estimate angular degree of freedom, we may set: ${\mathsf{\Omega }}_{\theta }\approx {\theta }^{\hspace{0.17em}3}/(8{\pi }^2)=4.1\times {10}^{-4}$, where $8{\pi }^2$ represents the maximum angular phase corresponding to a freely-rotating particle, and $\theta$ is the geometric mean of the three Euler angles accessible through trough particle rotation in the crystal.[59,63] We then extract $\theta \approx$18º, which is a physically acceptable[59] angular degree of freedom.

3.2. Thermodynamic Model for the Fluid Phase

In order to describe the thermodynamic behavior of protein solutions, an expression for $a_{\mathrm{R}}$ in Eq. 3 is needed. To enable LLPS, this expression must incorporate protein-protein attraction energy. According to our solubility results, HEPES is responsible for an increase in protein-protein attraction energy. This increase, however, should also cause a corresponding increase in LLPS temperature, in contrast with experimental results showing that HEPES causes a ≈5% decrease[52,53] in LLPS temperature. To further examine the effect of HEPES on protein-protein interactions, we examine previous measurements of lysozyme diffusion coefficient[52] as a function of protein concentration for the HEPES and REF models at 298 K. Assuming that the hydrodynamic factor[64] is the same for both systems, these diffusion data show that HEPES causes an increase of $1.2\pm 0.4$ in the second virial coefficient[43] ($B$ in ${\widehat{\pi }}_{\mathrm{HS}}/\varphi =1+B\hspace{0.17em}\varphi +\mathrm{...}$). This confirms that HEPES decreases protein-protein attraction energy in the fluid phase. It is in apparent conflict with HEPES being able to increase protein-protein attraction energy in in the crystalline phase.

We have previously mentioned that HEPES can stabilize lysozyme crystals by physical cross-linking interactions. This type of interaction can also occur in the fluid phase. However, their highly directional character makes them entropically unfavored, i.e., crosslinking can occur only if the relative orientation between two proteins is appropriate. In the crystalline phase, this entropy cost is marginal because protein orientation is essentially fixed by the lattice structure. In contrast, protein-protein attraction energy in the fluid phase can be weaker on average due to free rotation. In other words, it is the anisotropic character of protein-protein interactions that may explain the opposite effects of HEPES on lysozyme solubility and LLPS.

To theoretically derive LLPS boundaries, we need to discuss $a_{\mathrm{R}}$ in Eq. 3. There are several Monte-Carlo studies accurately describing the Helmholtz free energy of model fluid systems appropriate for protein solutions.[31,63,65,66,67,68,69] Nonetheless, it remains practically convenient to use thermodynamic-perturbation theories[44,58,59,70,71,72] that provide approximate analytical expression of the Helmholtz free energy. In our case, we shall consider a model that depends on just three parameters, describing attraction energy, range of interactions and degree of anisotropy. These are linked to the three main topological features of a dome-shaped LLPS boundary: the two critical point coordinates, $({\varphi }_{\mathrm{c}},{{\beta }_{\mathrm{c}}}^{-1})$ and the boundary width.

In this work, we shall consider the Barker-Henderson perturbation theory of square-well fluids (BHPT)[58,72] to describe the effect of HEPES on lysozyme phase behavior. Note that BHPT treats particle-particle interactions as isotropic. However, anisotropy may be incorporated in BHPT by introducing temperature-dependent energy parameter that depends on the fraction of protein surface engaging in protein-protein interactions.[63] Our revised BHPT model will be further discussed below. It is important to note that Wertheim perturbation theory of associating spheres (WPT)[44,59,70,73] has also been also applied to examine the phase behavior of protein solutions. WPT explicitly describes anisotropic interactions by considering a specified number of sites that can link two spherical particles. However, we chose to use BHPT instead of WPT as there is a noticeable difference in ${\varphi }_{\mathrm{c}}$ between WPT and computer simulation data[31,69] on spherical particles with same number of sites. Moreover, it also yields a relatively large difference in second virial coefficient between HEPES and REF systems. Nonetheless, application of WPT with the same number of parameters employed in our revised BHPT is discussed in Supplementary Material for completeness.

The residual Helmholtz free energy for the original BHPT model is written as

$$a_{\mathrm{R}}(\beta ,\varphi )=a_{\mathrm{mf}}+a_{\mathrm{R}2}$$

(8)

where $a_{\mathrm{mf}}(\beta ,\varphi )$ rigorously represents a mean-field first-order correction to the hard-sphere reference model. In the case of isotropic square-well potential, we have:[58]

$$a_{\mathrm{mf}}=-{\displaystyle \frac{1}{2}}\hspace{0.17em}{\nu }_{\mathrm{HS}}\hspace{0.17em}\varphi \hspace{0.17em}\beta \hspace{0.17em}{\epsilon }_{\mathrm{SW}}$$

(9)

where ${\epsilon }_{\mathrm{SW}}$ is the magnitude of the square-well attraction energy. Attraction occurs when the distance between the centers of two particles, $r$, is such that $\sigma \le r\le \lambda \sigma$ while no interactions occur if $r>\lambda \sigma$, where $\sigma$ is the particle diameter and $\lambda$ specifies the range of attractive interactions. In Eq. 9, $\hspace{0.17em}{\nu }_{\mathrm{HS}}$ is the average number of contacts that each particle makes within the range $1\le x\le \lambda$, evaluated using the pair distribution function of the reference hard-sphere fluid, $g_{\mathrm{HS}}(x,\varphi )$ with $x\equiv r/\sigma$. We specifically write:

$${\nu }_{\mathrm{HS}}(\varphi ,\lambda )=24\varphi {\displaystyle {\int }_1^{\lambda }{g}_{\mathrm{HS}}(x,\varphi )x^{2}dx}=8\hspace{0.17em}({\lambda }^3-1)\varphi \cdot g_{\mathrm{HS}}(1,{\varphi }^{\prime })$$

(10)

with $8\hspace{0.17em}({\lambda }^3-1)\varphi$ represents ${\nu }_{\mathrm{HS}}$ in the limit of $\varphi \to 0$ and $g_{\mathrm{HS}}(1,{\varphi }^{\prime })=(1-{\varphi }^{\prime }/2)/{(1-{\varphi }^{\prime })}^3$ is the Carnahan-Staling contact-value of $g_{\mathrm{HS}}(x,{\varphi }^{\prime })$ evaluated at the particle volume fraction, ${\varphi }^{\prime }(\varphi ,\lambda )$. The following Pade` approximant of ${\varphi }^{\prime }(\varphi ,\lambda )$ is available: [72,74,75]

$${\varphi }^{\prime }={\displaystyle \frac{c_1+c_2\hspace{0.17em}\varphi }{{(1+c_3\hspace{0.17em}\varphi )}^3}}\varphi$$

(11)

where $c_i={\displaystyle {\sum }_{j=1}^4{s}_{ij}{\lambda }^{-j}}$ with $i=1,2,3$ and the $3\times 4$ matrix of $s_{ij}$ coefficients is given by [-3.1649, 13.3501, -14.8057, 5.7029; 43.0042, -191.6623, 273.8968, -128.9334; 65.0419, -266.4627, 361.0431, -162.6996].

The second contribution in Eq. 8, $a_{\mathrm{R}2}$, takes into account higher order terms in BHPT and can be evaluated only approximately. The most successful approximation is the local compressibility approximation,[58,72] where $a_{\mathrm{R}2}$ is the second-order perturbation term of the Helmholtz free energy. The term, $a_{\mathrm{R}2}$, is approximately obtained by examining the radial profile of fluctuations of particle density around a central particle.[58] At a given distance from this particle, these fluctuations are locally described by the isothermal compressibility of the hard-sphere reference model. This leads to[58,72]

$$a_{\mathrm{R}2}=-{\displaystyle \frac{1}{2}}{\varphi }^2{\displaystyle \frac{d{\nu }_{\mathrm{HS}}}{d{\widehat{\pi }}_{\mathrm{HS}}}}{\displaystyle \frac{{\beta }^2\hspace{0.17em}{{\epsilon }_{\mathrm{SW}}}^2}{2}}$$

(12)

where $d{\nu }_{\mathrm{HS}}/d{\widehat{\pi }}_{\mathrm{HS}}=(d{\nu }_{\mathrm{HS}}/d\varphi )(d\varphi /d{\widehat{\pi }}_{\mathrm{HS}})$ and ${\beta }^2\hspace{0.17em}{{\epsilon }_{\mathrm{SW}}}^2/2$ represents the second-order term in the series expansion of the Mayer function, $e^{\beta {\epsilon }_{\mathrm{SW}}}-1$. We can numerically calculate $d{\nu }_{\mathrm{HS}}/d\varphi$ from Eq. 10 while we can use Carnahan-Starling equation of state[57,72] to determine that$d\varphi /d{\widehat{\pi }}_{\mathrm{HS}}={(1-\varphi )}^4/(1+4\varphi +4{\varphi }^2-4{\varphi }^3+{\varphi }^4)$. It is important to note that the BHPT model is applicable for moderately short interaction ranges; $\lambda \approx 1.2$ and higher.[75] This is appropriate for protein solution with ${\varphi }_{\mathrm{c}}\approx 0.22$ and less.[31] It should not be used when particle-particle interactions are very short ($\lambda <1.2$).

The range of interaction, $\lambda$, is the only parameter that affects the value of the critical volume fraction, ${\varphi }_{\mathrm{c}}$. In our application of the BHPT model, we should consider $\lambda$ as a fitting parameter that is just needed to reproduce experimental values of ${\varphi }_{\mathrm{c}}$. It is also important to note that the use of protein specific volume to convert concentrations into volume fractions is debatable as the diameter of a protein is not straightforwardly connected to protein specific volume. This adds uncertainty in the comparison between theoretical and experimental values of ${\varphi }_{\mathrm{c}}$. In some studies,[44,73] the protein diameter has also been proposed as an extra fitting parameter to be varied in order to enhance model accuracy.

It is known that the range of interactions specifies the value of ${\varphi }_{\mathrm{c}}$, with $\lambda \approx 1.3$ corresponding to ${\varphi }_{\mathrm{c}}\approx 0.20$.[31] Thus, we focus on the LLPS boundary extracted using Eqs. 8-12 with $\lambda =1.3$. For comparison, the critical point extracted from related Monte Carlo simulations[31] is also included. We can appreciate that the values of ${{\beta }_{\mathrm{c}}}^{-1}$ and ${\varphi }_{\mathrm{c}}$ calculated using BHTP are both ≈6% higher than those extracted from simulations. For completeness, we also included the boundary calculated after setting $a_{\mathrm{R}2}=$ 0. This exhibits relatively larger deviations from simulation results as expected.

It is important to note that Eq. 12 does not yield the correct expression of the second-virial coefficient, $B$, for the square-well fluid.[75] Specifically, it yields $B/4=1-({\lambda }^3-1)(\beta {\epsilon }_{\mathrm{SW}}+{\beta }^2{{\epsilon }_{\mathrm{SW}}}^2/2)]$ instead of $B/4=1-({\lambda }^3-1)(e^{\beta {\epsilon }_{\mathrm{SW}}}-1)]$. Since the second-virial coefficient is a fundamental thermodynamic property of fluids, it is important that the expression of $a_{\mathrm{R}2}$ reproduces the correct expression of $B$. We therefore replace Eq. 12 with

$$a_{\mathrm{R}2}=-{\displaystyle \frac{1}{2}}{\varphi }^2{\displaystyle \frac{d{\nu }_{\mathrm{HS}}}{d{\pi }_{\mathrm{HS}}}}(e^{\beta {\epsilon }_{\mathrm{SW}}}-1-\beta {\epsilon }_{\mathrm{SW}})$$

(13)

which is the same as Eq. 12 to second order in $\beta {\epsilon }_{\mathrm{SW}}$. As shown in Figure 4A, Eq. 13 significantly improves agreement with Monte Carlo Simulations. Specifically, the values of ${{\beta }_{\mathrm{c}}}^{-1}$ and ${\varphi }_{\mathrm{c}}$ are now found to be just ≈2% higher than those extracted from simulations. Hence, we will use Eq. 13 instead of Eq. 12 in our data analysis. In Figure 4B, we show the LLPS boundaries calculated at different values of $\lambda$. Here, we can see that ${\varphi }_{\mathrm{c}}$ increases as $\lambda$ decreases as expected from theory and Monte Carlo simulations.[31]

We now need to incorporate anisotropy in the BHPT model. One way to approximately describe anisotropy in the square-well fluid is to replace ${\epsilon }_{\mathrm{SW}}$ with a temperature-dependent effective energy, ${\epsilon }_{\mathrm{eff}}$, given by the following free-energy expression:[63]

$$\beta \hspace{0.17em}{\epsilon }_{\mathrm{eff}}/2=\hspace{0.17em}\ln [\alpha \hspace{0.17em}{e}^{\hspace{0.17em}\beta \hspace{0.17em}{\epsilon }_{\mathrm{SW}}/2}+(1-\alpha )]$$

(15)

where the factor “2” considers that a protein-protein interaction involves two proteins, $\alpha$ is a degeneracy factor representing the fraction of particle surface occupied by binding sites with energy ${\epsilon }_{\mathrm{SW}}$, and $1-\alpha$ represents the fraction of particle surface with zero binding energy. If $\alpha =1$, we recover the isotropic case: ${\epsilon }_{\mathrm{eff}}={\epsilon }_{\mathrm{SW}}$. As $\alpha$ decreases, the area covered by binding sites decreases and the anisotropic character of protein-protein interactions increases.

3.3. Effect of HEPES on Lysozyme LLPS Boundary

We shall now employ the BHTP model with $\lambda \hspace{0.17em}=\hspace{0.17em}1.3$ to describe the experimental LLPS boundaries for REF and HEPES systems. These are shown in the phase diagrams of Figure 5 together with the corresponding crystal solubility curves. To make the crystal model consistent with the BHPT model, we set the attraction energy parameter, ${\epsilon }_{\mathrm{SW}}$, for the fluid phase to be the same as ${\epsilon }_S$ extracted from solubility data. We are then left with identifying the values of $\alpha$that best fit the two sets of LLPS data. We determine that $\alpha =0.037$describes the LLPS of the REF system well. This value must be then decreased to $\alpha =0.027$ in order to accurately represent the lower LLPS temperatures of the HEPES system. Noteworthily, these values of $\alpha$ also describe the width of the LLPS boundaries fairly well. Specifically, experimental LLPS data show a decrease of ≈4% in LLPS temperature when protein volume fraction changes from ≈0.2 to 0.05. These variations are comparable with those shown for $\alpha =0.03$ in Figure 4D. In Figure 5, theoretical solubility curves are now obtained by calculating the protein chemical potential using the BHPT model. These theoretical curves, which are also shown in Figure 5, accurately describe solubility data justifying the initial use of $\widehat{\mu }\approx {\widehat{\mu }}_0+\ln \varphi$. Finally, we can also make predictions on the values of second virial coefficient at 298 K. Specifically, we use Eq. 15 to calculate that ${\epsilon }_{\mathrm{eff}}=$278 K for the REF system is higher than ${\epsilon }_{\mathrm{eff}}=$244 K for the HEPES system. If then use the square-well expression of second virial coefficient with ${\epsilon }_{\mathrm{eff}}$ replacing ${\epsilon }_{\mathrm{SW}}$:

$$B/4=1-({\lambda }^3-1)(e^{\hspace{0.17em}\beta \hspace{0.17em}{\epsilon }_{\mathrm{eff}}}-1)$$

(16)

we calculate $B=$-3.4 for the REF system and the higher value of $B=$-2.1 for the REF system. Their difference, +1.3, is essentially the same as that estimated from DLS data at 298 K ($+1.2\pm 0.4$).

4. Conclusions

We applied a colloid model to theoretically examine the opposite effects of HEPES on lysozyme crystal solubility and LLPS boundary. We specifically applied cell theory to solubility data to determine that HEPES increases protein-protein attraction energy by 2.7%. To explain the observed decrease in LLPS temperature, we considered that HEPES also enhances the anisotropic character of protein-protein interactions. We developed an analytical model based on BHPT to describe the LLPS boundaries using the same protein-protein attraction energy parameter extracted from solubility data (${\epsilon }_{\mathrm{SW}}={\epsilon }_S$). Only the parameter, $\alpha$, was decreased in order to increase anisotropy and successfully explain the effect of HEPES on LLPS temperature and change in second virial coefficient. Our work describes a useful analytical model for describing the multifaceted effects of additives on the phase diagram of protein solutions.