Simple Calculation of Phase Diagrams for Liquid-Liquid Phase Separation in Solutions of Two Macromolecular Solute Species.

A simple method is described for the calculation of two- and three-dimensional phase diagrams describing stability and coexistence curves or surfaces separating one- and two-phase regions in composition/temperature space of a solution containing solute species 1 and 2. The calculation requires a quantitative description of the intermolecular potentials of mean force acting between like (1-1 and 2-2) and unlike (1-2) species. Example calculations are carried out for solutions of species interacting via spherically symmetric square-well potentials as first-order models for protein-protein interaction. When the interaction between species 1 and 2 is more repulsive than those acting between like species, the two-phase region is characterized by an equilibrium between a phase enriched in 1 and depleted in 2 and a phase enriched in 2 and depleted in 1. When the interaction between species 1 and 2 is more attractive than those acting between like species, the two-phase region is characterized by an equilibrium between a phase enriched in both species and a phase depleted in both species. The latter example provides a first-order description of coacervate formation without postulating specific interactions between the two solute species.

Introduction

The study of liquid–liquid phase separation (LLPS) in solution dates back to classical investigations of the solution properties of synthetic and natural polymers.[1,2] Renewed interest in the subject has recently been engendered by the discovery of membraneless organelles within living cells that have been tentatively identified as immiscible liquid phases.[3−6] The membraneless organelles contain high concentrations of one or more proteins and/or nucleic acids and are thought to provide special microenvironments in which the rates and equilibria of critical biochemical reactions may be modulated or for sequestration of toxic substances.[5,7]

The present work concerns two classes of liquid–liquid phase transitions in solutions of two solute species, which we shall call species 1 and 2. In the first class, termed segregative transitions,[8] one phase is enriched (relative to the total composition) in species 1 and depleted (relative to the total composition) in species 2, which the second phase is enriched in species 2 and depleted in species 1. In the second class, termed associative transitions,[8] one phase is enriched in both solute species and the second phase is depleted in both solute species. The liquid phase that is enriched in both solutes is commonly referred to as a complex coacervate.[1]

Most theoretical analyses of the thermodynamics of LLPS in solutions of two macromolecular solute species follow the general approach pioneered by Flory and colleagues,[2] according to which two immiscible solution phases form when the free energy of mixing the two solute species is positive. The free energy of mixing is decomposed into enthalpic (or energetic) and entropic contributions, and models for the composition dependence of each contribution are specified (see for example refs (9) and (10)). A second general approach was proposed by Edmond and Ogston,[11] according to which the composition dependence of the chemical potentials of each of the two species is calculated or measured, and the composition of phase boundaries obtained by satisfying conditions for thermodynamic stability and equilibrium. Analyses of this type have been employed in studies of the properties of solutions of mixed crystallins.[12−14]

The major difference between the two approaches is that the latter does not require specification or quantitation of individual enthalpic and entropic contributions to the total free energy of the solution. It is therefore more general in principle, since it does not distinguish between different types of biological macromolecules (proteins, nucleic acids, and polysaccharides) mixtures of which would require qualitatively different models for the enthalpy and entropy of mixing. Moreover, as will be elaborated below, the composition dependence of the chemical potentials of each solute species in a mixture may be experimentally measured in a straightforward fashion via analysis of the composition dependence of colligative properties of solutions of individual species and their mixtures. For this reason, in the present work we shall adopt and extend the free-energy-based approach of Edmond and Ogston,[11] and will demonstrate that it is capable of exhibiting equilibrium behavior that is characteristic of both segregative and associative phase transitions.

In the following section, the calculation of composition-dependent chemical potentials is described. Next the calculation of composition- and temperature-dependent stability boundaries (spinodals) and equilibrium phase boundaries (binodals) and determination of solution composition within the two-phase region is described. Then, results of example calculations exhibiting properties of both segregative and associative LLPS are presented. Finally, the significance of these results is discussed, and comparisons are drawn with other treatments of LLPS.

Methods

Specification of Solute Chemical Potentials

The chemical potential, or free energy increment, of solute species i in a solution of an arbitrary number of solute species may be generally written aswhere μ0 denotes the standard state chemical potential, R is the molar gas constant, T is the absolute temperature, c is the molar concentration of species i, γ is the thermodynamic activity coefficient of species i, and {c} is the concentrations of all solute species. ln γ is a direct measure of the equilibrium average free energy of interaction per mole between a solute molecule of species i and all the other solute molecules in the solution. As will be seen below, this quantity and its derivatives with respect to solute concentrations will be utilized to directly calculate LLPS phase boundaries.

The first step in simplifying our approach to the analysis of phase equilibria in solutions is recognition that the mathematical description of the equilibrium thermodynamic properties of solutions is isomorphic to that of a fluid of solute particles,[15] provided that the potentials of direct interaction (i.e., through vacuum) between the solute particles in the fluid are replaced by a potential of mean force acting between solute molecules in solution, defined as follows for the interaction between two molecules:The first term on the right side of eq denotes the free energy of the solution when the centers of molecules i and j are separated by distance r with a mutual orientation denoted by the generalized variable ω, and the second term denotes the free energy of the solution when the molecules are separated by a sufficiently large distance such that the free energy of the solution is no longer a function of r. The potential of mean force acting between solute molecules i and j in solution thus incorporates not only direct interaction between the two solute molecules but also contributions to the free energy of interaction due to solvent and all solute species except molecules i and j. According to the McMillan–Mayer theory of solutions,[15] the thermodynamic activity coefficients of each of two solvent species in solution may be expanded in powers of solute concentration according towhere the two-body interaction coefficients B are defined as functions of the potential of mean force acting between molecules of species i and j, and the three-body interaction coefficients B are defined as functions of the potential of mean force acting between molecules of species i, j, and k. When truncated after three-body terms, eqs and 4 are valid only over a restricted range of concentrations, the upper limit of which depends upon the strength of the intermolecular interactions under a particular set of experimental conditions.

There are three unique two-body interaction coefficients and four unique three-body interaction coefficients. These seven interaction coefficients are the only input required in order to calculate phase diagrams, provided that the concentration is sufficiently low that truncation of eqs and 4 after three-body terms is valid. They may be evaluated experimentally via measurement of the composition dependence of any of three thermodynamically determined colligative properties of solution mixtures–osmotic pressure,[16] sedimentation equilibrium,[17] and light scattering.[18,19] They may also be evaluated computationally given suitable relations for the potential of mean force acting between two and three solute molecules in solution as a function of the relative distances and orientations of the interacting molecules.[13,14,20] A simplified version of this last option is presented below.

Specification of Potential of Mean Force Acting between Solute Molecules

In order to demonstrate our simplified approach to the calculation of LLPS phase diagrams, we shall evaluate the B and B, and via eqs and 4, the chemical potentials of each of two solute species in a mixture interacting via square-well potentials of mean force using an analytical formalism due to Kihara.[21,22] Relations derived by Kihara, specifying the interaction-dependent values of B and B, are presented in the appendix to ref (23). The three two-body and four three-body interaction coefficients so calculated have been found to accurately describe the activity coefficients of each solute species as a function of composition in mixtures of model square-well interacting solutes at total volume fractions of up to 0.2, as calculated via Monte Carlo simulations.[23] It is evident that the SW potential is only a crude approximation to any real potential of mean force acting between real biological macromolecules. However, it does incorporate the two major features of any real potential of mean force: a “hard” steric repulsion and a “soft” longer-ranged attraction or repulsion. Moreover, it will be shown subsequently that results obtained with this very simple model can reproduce qualitatively the phase separation behavior exhibited by real protein and polymer solutions.

For the purpose of demonstration, we define the steric radius of an effectively spherical molecule to be equal to the radius of a hard sphere with the same mass and density as the actual molecule.where M denotes the molar mass of species i, v̅ is the partial specific volume, and NA is Avogadro’s number. The square-well potential of interaction between molecules of species i and j is then defined as follows:where r denotes the distance between the centers of molecules of species i and j, ε is the depth (or height) of the square well, k is the Boltzmann constant, and L is a scaling parameter relating the range of the square-well interaction to the sizes of the two interacting molecules. For convenience we shall subsequently introduce the variable ε* ≡ ε/kT to denote the value of ε in units of the thermal energy kT. This potential is shown schematically in Figure . Effective interactions between two species of solute molecules in solution are thus defined by the parameters r1, r2, ε11*, ε22*, ε12*, L11, L22, and L12. In the calculations to follow we shall assume for the purposes of demonstration that the partial specific volumes of the interacting macromolecules are equal to 0.73 cm3/g, an average value for globular proteins,[24] so the values of r1 and r2 are determined by the corresponding molar masses.

Figure 1

Plot of square-well potential of mean force as a function of the center-to-center distance between two interacting spherical molecules of species i and j.

Determination of Phase Boundaries

Stability Boundaries (Spinodals)

Compositions lying on the spinodal are characterized by a free energy surface that is concave down in at least one direction in composition space, so any small fluctuation in composition will lead to a rapid phase separation termed spinodal decomposition. According to the fluctuation theory of light scattering in a solution of two scattering solute species,[18,25] at spinodal compositions the scattering intensity will diverge, leading to the conditionEvaluation of the composition-dependent activity coefficients of species 1 and 2 and their partial derivatives via eqs and 4 together with the numeric solution of eq permits calculation of the spinodal curves.

Equilibrium Phase Boundaries (Binodals)

Compositions lying on an equilibrium boundary must satisfy the condition that the chemical potential, or free energy increment, of each solute species must be equal in both phases.where the superscript indicates the phase. Since the standard state chemical potential of each species is independent of composition and therefore identical in both phases, eqs and 9 are equivalent towhere Δμ ≡ (μ – μ0)/RT = ln c ln γ(T, {c}). Solution of these equations leads to a set of pairs of equilibrium compositions (c1I, c2I) and (c1II, c2II), which may be connected by straight tie-lines. Any total composition (c1, c2, T) within the two-phase region will lie on a tie-line connecting one pair of these equilibrium compositions, and the volume fraction of phase II will be given byThe numeric solution of eqs and 11 for the entire set of equilibrium concentrations is prohibitively computationally intensive in the general case. For example, if one wished to explore a range of concentrations in w1 and w2 with a resolution of 100 points in each concentration range, one would have to evaluate chemical potentials for 1004 combinations of concentrations. In the present work, we restrict calculation of binodals to the symmetric case: M1= M2, L11= L22, and ε11* = ε22*, permitting calculational shortcuts as described in the Supporting Information.

Simulating Temperature Dependence of Phase Transitions

We define a reference temperature T0, a reference value of ε* = ε/kT0, and the relative temperature Trel = T/T0. In principle, one should allow for the temperature-dependence of intermolecular interactions. Without specifying the nature of these interactions, we may writewhere the value of α indicates whether the strength of the underlying interaction increases (α > 0) or decreases (α < 0) with temperature. In order to simulate the temperature dependence of phase transitions, values of α and ε* are specified, corresponding to a relative temperature of unity. Then binodal compositions are calculated as described above for each of an array of increasing values of Trel. If α < 1, then the absolute value of ε*(Trel) diminishes with increasing Trel, such that at some maximum value of T that we designate the critical value TUC, or upper consolute temperature, the phase boundary will vanish, and the solution will remain a single phase at all temperatures greater than Tcrit. If, however, α > 1, then the absolute value of ε*(Trel) will increase with increasing Trel, in which case a solution existing as a single phase at Trel = 1 may form two phases above a critical temperature T, referred to as a lower consolute temperature (see below).

Numerical Methods

Numeric solution of eqs , 4, 7, 10, and 11 was performed using user-written scripts and functions in MATLAB (Mathworks, Natick MA) that are available upon request from the author. The algorithms employed are described in the Supporting Information.

Results

Segregative Phase Transition

Figure shows a phase diagram obtained from a calculation carried out at a single temperature with the parameter values specified in the caption. Note that the value of ε12* is more positive than either ε11*or ε22*, indicating that the interaction between molecules of unlike species is more repulsive than the interaction between molecules of either like species. The blue curve indicates the spinodal, and the red curve indicates the binodal. The dashed tie-line connects the compositions of two equilibrium phases: (w1I = 224,w2I = 33.5) and (w1II = 33.5,w2II = 224). Symmetry dictates that the chemical potentials of each solute is the same in both phases. Qualitatively similar results were obtained from other simulations in which the interaction between unlike species is more repulsive than the interactions between like species.

Figure 2

Phase diagram at fixed temperature for symmetric segregative interactions between two solute species. Results of calculations carried out for M1 = M2 = 70 000, L11 = L22 = L12 = 1.25, ε11* = −0.5, ε22*= −0.5, and ε12* = 2. The spinodal is plotted in blue and the binodal in red. The dashed tie-line connects two phases with equilibrium compositions (224, 33.5) and (33.5, 224).

An experimental measurement of phase equilibria in a solution mixture of two polymers[11] yields a phase diagram, shown in Figure , that resembles qualitatively the calculated binodal in Figure .

Figure 3

Binodal curve separating one-phase region (lower left) from two-phase region (upper right) in a solution mixture of PEG-6000 (species 1) and Dextran 19.7 (species 2). Data from Edmond and Ogston.[11]

Figure shows the results of a calculation carried out at multiple temperatures as described above, with parameter values specified in the figure caption. The left panel is a projection of the three-dimensional phase diagram onto the composition axis. The blue curves are binodals calculated at various relative temperatures plotted as functions of solution composition. Since α is set equal to 0, the absolute values of ε*(Trel) diminish with increasing Trel. The left panel shows that with increasing temperature, binodals trend toward the upper right-hand corner. The red curve defines the intersection between the binodal surface and a vertical plane connecting compositions denoted by (w1, w2) of (0, 240) and (240, 0), indicating the critical temperature for a solution of the corresponding composition. At temperatures exceeding this value the solution exists as a single phase, and a temperatures below this value, at equilibrium the solution exists as a mixture of two phases. As an example, the horizontal tie-line plotted at Trel = 1.18 connects two phases with compositions (59.9, 108.5) and (180.5, 59.9). The black dot indicates the value of the critical temperature at a solution composition of (120, 120).

Figure 4

Temperature-dependent phase diagram for symmetric segregative interactions between two solute species. Results of calculations carried out for M1 = M2 = 70 000, L11 = L22 = L12 = 1.25, ε11,0* = ε22,0*= −0.5, ε12,0* = 2, and α = 0. Left panel: projection of the three-dimensional binodal surface onto the composition axis. Right panel: binodals calculated for various relative temperatures plotted as a function of composition. The red curve represents the intersection of the binodal surface with a vertical plane extending between (0, 240) and (240, 0), the meaning of which is described in the text.

An experimental measurement of phase equilibria in a solution mixture of two polymers carried out at multiple temperatures[26] yields a three-dimensional phase diagram, shown in Figure , that qualitatively resembles the calculated diagram in Figure .

Figure 5

Composition–temperature phase diagram for a solution mixture of soluble gelatin (weight-average molar mass = 170 kg) and dextran (weight-average molar mass = 282 kg). Data from Edelman et al.[26] Reproduced (with modification) with permission from ref (26). Copyright 2001 American Chemical Society.

Associative Phase Transition

Figure shows a phase diagram calculated at a single temperature with the parameter values shown in the figure caption. The value of ε12* is more negative than either ε11*or ε22*, indicating that the interaction between molecules of unlike species is more attractive than the interaction between molecules of either like species. The left panel shows closed binodal (red) and spinodal (blue) curves separating an exterior one-phase region and an interior two-phase region. The dashed spinodal curve is calculated via numerical solution of eq , as described in the Supporting Information. The solid spinodal and the binodal curves are calculated using an approximate analysis of the composition dependence of osmotic pressure, as described in Supporting Information. Two sample tie-lines are plotted, each of which is connects points indicating phase compositions that are depleted in both solute species and enriched in both solute species, relative to any total composition lying along the tie-line. For example, the uppermost tie-line plotted in red connects compositions of (21.1, 32.2) and (40.5, 55.4). In the right panel, the (relative) chemical potentials of each of the two solute species are plotted as functions of binodal composition, and it may be seen how the two compositions indicated by the termini of each tie-line in the left-hand panel connect points of equal chemical potential of each species in both phases, satisfying eqs and 12. Qualitatively similar results were obtained from other simulations in which the interaction between unlike species is more attractive than the interactions between like species.

Figure 6

Results of calculations carried out for M1 = M2 = 70 000, L11 = L22 = L12 = 1.25, ε11* = ε22* = −0.8, and ε12* = −2.5. Left panel: phase diagram at a fixed temperature showing calculated spinodal (blue) and binodal (red), with two sample tie-lines. Right panel: relative chemical potentials of each solute species plotted as a function of composition along the binodal curve. Points are plotted for termini of each of the tie-lines shown in the left panel.

An experimental measurement of phase equilibria in a solution mixture of two polymers[1] yields a phase diagram, shown in Figure , that resembles qualitatively the calculated diagram in the left panel of Figure .

Figure 7

Tie-lines connecting compositions of equilibrium phases in a solution mixture of gum arabic (species 1) and gelatin (species 2). Symbols: data from Figure 6 of Bungenberg de Jong.[1] Blue dashed curve: smooth extrapolation of the data representing a possible binodal curve connecting the data points.

Figure shows a temperature-dependent phase diagram calculated as described above, with parameter values specified in the figure caption. In the left panel the red and blue curves are binodals and spinodals, respectively, calculated for each of a series of values of Trel, plotted as functions of composition. As the temperature increases, the size of the two-phase region decreases and ultimately vanishes when the temperature exceeds Tcrit, the value of T above which the solution remains one-phase at all compositions. A symmetry axis defined by w1 = w2 = 0.5 wtot is drawn in black, connecting the compositions of two phases coexisting at Trel = 1.04. The combination of constant-temperature binodals defines a three-dimensional region of two phases in composition/temperature space. In the right panel, points are plotted to indicate the two intersections of the binodal and two intersections of the spinodal at each temperature along the symmetry axis. These points define composition–temperature binodal and spinodal curves in the plane of the equal-concentration slice.

Figure 8

Simulated temperature dependent phase diagram. Blue curves are spinodals and red curves are binodals, calculated as described in the text, with M1 = M2 = 70 000, L11 = L22 = L12 = 1.25, α = 0, ε11,0* = ε22,0* = −0.8, and ε12,0*= −2.5. The dotted black line in the left panel is the axis of symmetry (w1 = w2), and the marked points on the binodal and spinodal curves represent the values of the binodal and spinodal compositions along the axis of symmetry. These points are replotted as a function of total composition in the right panel.

While we have so far been unable to find a published experimental measurement of associative LLPS at multiple temperatures with which to compare our calculated results, Figure shows experimental results obtained at multiple salt concentrations. Increasing salt concentration, like increasing temperature, will diminish the strength of electrostatic intermolecular interactions, and one would expect the effect of increasing salt concentration to qualitatively resemble the effect of increasing temperature.

Figure 9

Compositions of phases in equilibrium in mixtures of gum arabic (species 1) and gelatin (species 2) at various concentrations of CaCl. Symbols: data from Figure 7 of Bungenberg de Jong.[1] Black dotted line indicates axis of equal w/v concentrations of the two polymers and connects compositions of equilibrium phases at 4% salt. Dashed red curves are smooth extrapolations of the data representing possible binodals. Right panel is a plot of the points at all salt concentrations along the equal w/v composition axis. Blue curve in right panel is a cubic polynomial fit to the points, plotted to guide the eye.

Discussion

The resemblance between the calculated results shown in the previous section and the experimental results with which they are compared is qualitative rather than quantitative. For several reasons, this is to be expected: (1) The square-well model for a potential of mean force is clearly simplistic and used only to demonstrate application of eqs ,4, and 7–9 to the calculation of spinodal and binodal curves. (2) Equation , used to simulate a temperature-dependent phase diagram respectively, is heuristic and meant only to capture qualitative effects. (3) The experimental results cited for comparison with calculations were obtained using mixtures of polymers that are quite heterogeneous rather than individual species as specified in the model. Nevertheless, the correspondence between calculated and experimental phase diagrams is clear. When interaction between unlike species is less favorable than that between like species, the predicted tie-lines connect phases enriched in one species and depleted in the other, characteristic of segregative LLPS, and when interaction between unlike species is more favorable than that between like species, the predicted tie-lines connect phases enriched in both species with phases depleted in both species, characteristic of associative LLPS or coacervate formation.

The phase diagrams shown in Figures and 5 exhibit upper consolute temperatures, indicating that a decrease in temperature favors phase separation. Such behavior is often observed,[8,10] but systems exhibiting lower consolute temperatures exist as well,[27] where phase separation is induced upon increasing temperature. The appearance of phase separation with increasing temperature indicates that intermolecular interactions are strengthened with increasing temperature. Such systems may be qualitatively simulated by the method introduced here, provided that calculations are performed using eq with a value of α exceeding 1, as shown in Supporting Information.

The phenomenon of LLPS has been observed in mixtures of two polymers,[1,8] mixtures of a polymer and a protein,[28−30] and in mixtures of two proteins.[14,31−33] The thermodynamic approach to calculation of LLPS in solutions of two solutes pioneered by Edmond and Ogston[11] and generalized here may be applied to the analysis of all of these phenomena, whereas the classical approach employing separate models for the energy (or enthalpy) and entropy of mixing of two solute species must be tailored specifically for each mixture. For example, the frequent utilization of Flory–Huggins mean-field theory[2] to estimate the entropy of mixing of two solutes clearly does not apply when one or both of the solutes are not polymers, and models proposed for interaction of different types of macromolecules clearly vary qualitatively with the nature of the macromolecules. In contrast, the values of the interaction coefficients B and B appearing in eqs and 4 comprise the totality of information required to perform the calculations presented here. The combination of eqs , 4, and 7 informs us that a spinodal boundary will occur at all solution compositions (c1, c2) satisfying the model-independent relationwhereAs pointed out above, the required values of the B and B may be obtained by means of model functions for the potentials of mean force, as in the present work and in refs (12−14) by molecular dynamics or Monte Carlo computer simulations,[34] or from experimental measurement of the composition dependence of colligative properties such as static light scattering,[18,19] sedimentation equilibrium,[17] or osmotic pressure.[16,35] While the methods of calculating binodals presented here is limited to symmetrically interacting systems, the calculation of spinodals via solution of eq is applicable to asymmetrically as well as symmetrically interacting systems, as shown in Figure . If an experimentally measured set of measured values of B and B results in a predicted spinodal, then one may be confident that the solution will exhibit a liquid–liquid phase separation when the total composition falls within the predicted spinodal boundary. In contrast, solutions with compositions falling between the binodal and spinodal may or may not separate into two phases, depending upon whether the solution has or has not achieved thermodynamic equilibrium.

Figure 10

Spinodal boundaries calculated for symmetrically and asymmetrically interacting solute mixtures by numerical solution of eq . Black spinodal is calculated for a symmetrically interacting mixture with the parameters given in the caption of Figure . The blue and red spinodals are calculated for asymmetrically interacting solute mixtures with the following interaction potentials. Blue: L as above. Black: ε11* = −0.7, ε22* = −0.8, and ε12* = −2.5. Red: L11 = 1.25, L22 = 1.28, L12 = 1.25, ε11* = −0.7, ε22* = −0.8, and ε12* = −2.5. Black dotted line is the symmetry axis.

Although the analysis of LLPS from the standpoint of chemical potential was pioneered by Edmond and Ogston (EO),[11] the work presented here is considerably more general than the original: (1) The present work takes into account three-body solute–solute interactions, whereas EO took into account only two-body interactions, and thus is applicable at higher solute concentrations. (2) Unlike EO, the present work provides a protocol for calculation of spinodals. (3) The present work provides methods for calculating phase-diagrams for associative LLPS, including both spinodals (in general) and binodals (in special cases), whereas EO did not treat associative LLPS at all.

At the phenomenological level considered here, solvent, small-molecule solutes such as salts, and temperature are treated implicitly rather than explicitly. However, in order to arrive at a quantitative mechanistic description of the effects of solvent, small molecules, or temperature on LLPS in a solution containing two specific macromolecular solutes, a quantitative analysis of the dependence of the potentials of mean force acting between molecules of like and unlike macromolecular solutes upon each of these factors will ultimately have to be incorporated. Such information may be obtained experimentally through measurements of the effect of varying temperature or salinity upon the solution properties listed above or via sufficiently detailed computer simulation.

Finally, we emphasize that the “simple” calculation of phase diagrams presented here is rapid. Using a generic desktop personal computer running MATLAB (Mathworks, Natick, MA), together with scripts and functions incorporating algorithms described in the Supporting Information, a phase diagram such as that shown in Figure may be generated in a few seconds and that shown in Figure within a minute. This rapidity allows the investigator to explore a large parameter space including variations in molecular size, and the strength and range of both attractive and repulsive interactions between like and unlike molecules in solution mixtures over a wide range of composition and temperature.