Phase behavior in multicomponent mixtures

In this article, we study the phase behavior of two polydisperse hydrocolloids: dextran and polyethylene oxide. We combine the data on the experimental osmometric virial coefficients of the pure components with the experimental critical point of their aqueous mixture and the size distribution of each component from a previously published study in order to predict the phase boundary, spinodal, and fractionation upon demixing of the polydisperse mixture. We compare the results of our calculation to the experimental phase diagram. Our method reveals a better correspondence with the experimental binary phase behavior than modeling each component as monodisperse. The polydispersity of the hydrocolloids causes the phase separation boundary to shift to lower concentrations and the miscibility region to decrease and change its shape from a rotated U-shape to a W-shape.

1 Introduction

A solution may gel when assemblies of molecules in that solution span the solution. The properties of the assemblies and their mutual interactions determine the concentration range above which the gel forms and the properties of the gel. For example, in the case of fibrillar assemblies from protein-derived peptides, gelation occurs already at low concentrations of the building blocks of the fibril (van der Linden, 2012). Notably, such fibrils exhibit interesting collective properties at higher concentrations, such as their arrangement in an anisotropic manner like that of (nematic) liquid crystals (Bagnani et al., 2019).

The molecular assembly may reflect an equilibrium state, but, for most practically relevant systems, in particular for gelling systems, it more often reflects a non-equilibrium state. The structure of such a state depends on the specific spatial–temporal path through equilibrium and unstable regions in the phase diagram that are followed during the preparation of the system and on the rate of change of the thermodynamic conditions. As such, the equilibrium phase diagram plays an important role in predicting the characteristics of the structures within and the concomitant properties of the gel systems. This importance of phase diagrams holds for systems of any number of components. We will discuss a few examples of the above in simple systems in this Introduction section to set the scene, before addressing phase behavior in multicomponent mixtures in more detail in the next section.

As a first example, we choose a system composed of one component in water. We choose gelatin, which forms random coils at high temperatures and is a liquid, while below a critical temperature, its molecules form triple helices and its system shows gelation at sufficiently large concentrations. The work of Chatellier et al. (1985) demonstrates that a critical concentration of helices exists, below which no gel forms. Joly-Duhamel et al. (2002a, 2002b), in addition, have clearly shown that the liquid–gel transition is, in fact, a percolation transition, with the concentration of triple helices being the control parameter. We have shown that the critical concentration of helices, below which the system gels, can be deduced from the structural characteristics of the triple helices, i.e., their persistence length and thickness (van der Linden and Parker, 2005). The helices can, in fact, be viewed as fibrils. Furthermore, the value of gel elasticity at low concentrations of helices can be solely deduced from entropy arguments (van der Linden, 2012), and this also holds for the elasticity value at higher concentrations (apart from the fact that one requires one adjustable fitting parameter which is predicted to be of order 1, as confirmed by fitting) (van der Linden, 2012). The elasticity over time is determined by the final temperature of the system and the rate at which this occurs (Normand et al., 2000; Normand and Parker, 2003). We clearly have a non-equilibrium gel state. Despite this non-equilibrium state, one can always predict elasticity solely on the basis of the helix concentration, which is a separately accessible observation of Joly-Duhamel et al. (2002a, 2002b).

For another case of fibrillar assembly under non-equilibrium conditions, Nguyen and Vaikuntanathan (2016) simulated the fibrilization of two components A and B for a given initial excess chemical potential of each type of monomer. This situation can be experimentally realized, for example, in a supersaturated solution of both types of monomers. The excess chemical potentials were assumed to be equal for both types. The interactions between the likewise species are also assumed to be equal. The fibrils were found to consist of blocks composed of A, connected to blocks composed of B. The less favorable the unlike particle interaction becomes, the longer the blocks become. If one now increases the initial excess chemical potential of the two types of monomers, the length of the blocks becomes smaller. As such, the structural diversity within a fibril increases upon increasing the initial excess chemical potential of the monomers, i.e., upon increasing the distance to the transition point. This in turn will affect the interactions and flexibility of the fibril and thus the properties of the gel it forms.

Another example of a “one-component” system, for which the equilibrium phase behavior is important in understanding the non-equilibrium gel state, consists of oil droplets in water, for which the interaction can be tuned by the temperature (Poulin et al., 1999). Above a sufficiently large enough attraction between the droplets, the authors find a transition, known as spinodal decomposition. This occurs in the region of the phase diagram that is referred to as the instability region. One finds a bi-continuous structure, with one of its parts being mostly droplets that are packed with a fractal dimension, d_f. This d_f can get as low as 1.7 at sufficiently high attraction, forming a gel state, while for lower attractions, one finds a d_f of approximately 3, indicating a fluid state, resulting in the end in a two-phase (de-mixed) system.

In order to better understand gelled systems made from multicomponent mixtures, the above urges us to investigate the phase diagrams of multicomponent mixtures. Hereto, we currently report on the recent theoretical and experimental results obtained to explain the experimental phase diagram for a multicomponent system in water containing differently sized coil-like polymers of dextran (D) and polyethylene oxide (PEO). Mixtures of such polydisperse bio-polymeric ingredients are ubiquitous and therefore bear relevance to many practical situations. An important example is the understanding of functionality in food formulations of a variety of multicomponent mixtures that consist of minimally refined plant-based ingredients. This understanding facilitates the adaptation to use different ingredient sources, thereby supporting the development of a more sustainable food supply.

Aqueous mixtures of polymers such as polyethylene oxide and dextran form liquid–liquid two-phase systems at certain concentrations. These systems are often used as a model system for the phase behavior of macromolecules in solution because they show a clear macroscopic phase separation (Kang and Sandler, 1988; Edelman et al., 2003; Dewi et al., 2020). They also have practical applications, as they are often used to aid in the partitioning of biological materials such as proteins and cell materials (Johansson et al., 1998; Johansson and Walter, 1999).

For the prediction of their phase behavior, the polymers are often considered monodisperse. Experimental work by Edelman et al. (2003), however, has shown that upon demixing, considerable fractionation in the molecular weight occurs for both PEO and dextran. Not only is there an effect of the polydispersity of each component on fractionation but also the concentration in the parent phase plays a role. At lower total polymer concentrations in the parent phase, the depleted colloids have a broader distribution, i.e., there is less pronounced fractionation. The changes in the distribution are mainly prevalent in the depleted phase. The average molar mass in the enriched phase does not change considerably. The changes in the amount of fractionation, depending on the concentration along a dilution line, were also reported by Zhao et al. (2016). Furthermore, Gaube et al. (1993) showed that polydispersity plays a role in phase behavior and phase composition. They compared mixtures of PEO and dextran with various molecular sizes and found that for each polymer, the short-chain molecules preferentially partition to the phase enriched in the other polymer.

There has been some effort in incorporating the polydispersity of these polymers in the prediction of their phase behavior. Using a universal quasi-chemical (UNIQUAC) model, Kang and Sandler (1988) incorporated polydispersity in their prediction and found that polydispersity of the polymers enlarged the two-phase region considerably near the critical point and resulted in smaller miscibility regions far from the critical point. They also found significant fractionation, and the difference in the average molecular weight of the components in each phase increased with larger polydispersity.

Most often, when the phase behavior of a binary mixture is studied experimentally, one or more dilution lines are used to obtain the concentration of each component in each phase. These concentrations are then used to construct the binodal (Albertsson, 1970). However, this approach does not shed light on polydispersity nor does it give an insight into the demixing in the metastable region, where a system de-mixes into multiple phases, but that is outside of the unstable (spinodal) region. It is also difficult to quantify the impact of polydispersity and its corresponding distribution on phase behavior (Dos Santos et al., 2004). In order to get more insights into the transition between one and two phases, Larsson and Mattiasson (1988) determined the experimental phase boundary for polydisperse PEO and dextran. They found a significant broadening of the phase boundary compared to the approach of Albertsson (1970). The broadening increased with increasing polydispersity at the depleted side of the said polydisperse component (Larsson and Mattiasson, 1992).

In a previous numerical work, we predicted the phase behavior of a polydisperse binary mixture of hard spheres in solution using the virial coefficient approach (Sturtewagen and van der Linden, 2021; Sturtewagen and van der Linden, 2022), based on the theory of McMillan and Mayer (1945). Polydispersity was incorporated into the system by means of sub-components. The spheres had an asymmetric size ratio to aid in demixing (Dijkstra et al., 1999). Polydispersity of the largest component caused significant changes to the phase diagram. With an increase in polydispersity, the critical point shifted to higher concentrations, while also at the same time, the miscible region decreased. Not only polydispersity plays a role but also the type of distribution has an influence. We found that the largest components of the distribution impacted the phase diagram the most. De-mixed mixtures also showed significant fractionation. The smallest sub-components of the large polydisperse components favored the top phase that was enriched in the small monodisperse component.

In recent works, the experimental phase diagram for the macromolecules of PEO and dextran was theoretically reconstructed using osmotic virial coefficients as obtained from the fitting experimental data from an earlier published work (Dewi et al., 2020) and the experimentally obtained critical point of the binary mixture (Dewi et al., 2020). However, some aspects of the experimental phase diagram remained unexplained by the theoretically constructed one. Mainly, at low concentrations of dextran, the experimental system showed thermodynamic incompatibility and demixing into two phases, while the theory predicted a homogeneous mixture for these concentrations. We hypothesize that this discrepancy is due to the polydispersity of the macromolecules.

In this article, we take the polydispersity of both these components into account (by subdividing each polymer into sub-components of different sizes) when predicting the phase diagram based on the experimental data from Dewi et al. (2020). We study the position of the spinodal and binodal. We model the interactions between the different polydisperse sub-components using the second virial coefficient, where we assume that the polydisperse sub-components act as non-additive hard spheres. We start with describing the underlying theoretical background in Section 2. In Section 3, we describe the materials and methods. In Section 4, we discuss the results and compare it with that of other works in the area of multicomponent mixtures.

2 Theoretical background

To describe phase separation in a system containing particles of different diameters and of different types, we choose a framework that accounts for terms up to quadratic terms in concentration, i.e., including up to the second-order virial coefficients that are experimentally or theoretically accessible. We start from the Helmholtz free energy $F$ for an N-component mixture in a common solvent:

$\frac{F}{RTV}={\displaystyle \sum _{i=1}^N}{c}_i\ln \left(c_i\right)+{\displaystyle \sum _{i,j=1}^N}{B}_{ij}{c}_i{c}_j,(1)$

where $B_{ij}\equiv B_{ji}$, $R$ represents the molar gas constant, $T$ represents the temperature, and $c_i$ represents the molar concentration of component i, with $c_i=n_i$/V, where $n_i$ denotes the number of moles of component i in volume $V$, and $B_{ij}$ represents the second virial coefficients, reflecting the interactions between components i and j. The osmotic pressure $\mathrm{\Pi }$ of a mixture can be written as

$\frac{\mathrm{\Pi }}{RT}=-\frac{1}{RT}{\left(\frac{\partial F}{\partial V}\right)}_{T,n_i}={\displaystyle \sum _{i=1}^N}{c}_i+{\displaystyle \sum _{i,j=1}^N}{B}_{ij}{c}_i{c}_j.(2)$

In our case, the interactions between the particles are assumed to follow the non-additive hard sphere model (Sturtewagen and van der Linden, 2022), where the non-additivity refers to either an attractive or repulsive extra interaction on top of the excluded volume interaction. For such non-additive hard spheres, the virial coefficients for a system of two compounds, $i$ and $j$, consisting of spheres of two different diameters, are given by

$B_{ij}=\frac{2}{3}\pi ·{\left[\frac{{\sigma }_i+{\sigma }_j}{2}·\left(1+{∆}_{ij}\right)\right]}^3,(3)$

where ${\sigma }_i$ refers to the diameter of the sphere of the compound i, and ${∆}_{ij}$ reflects the interaction contribution that is different from the excluded volume interaction. The term ${∆}_{ij}$ can be negative or positive in the non-additive regime, whereas it is zero while referring to a pure excluded volume interaction, with the latter referring to the so-called additive hard sphere model. The chemical potential µ_i for component i can be obtained from

$\frac{{\mu }_i}{RT}={\frac{1}{RT}\left(\frac{\partial F}{\partial n_i}\right)}_{T,V,n_{i\ne j}}=\frac{{{\mu }_i}^0}{RT}+\mathit{ln}{c}_i+2{\displaystyle \sum _{j=1}^N}{B}_{ij}{c}_j,(4)$

where ${{\mu }_i}^0$ is the standard chemical potential of compound $i$. Suppose we search for conditions for one or two macroscopic phases, there are three important aspects: the occurrence of instability, the critical point, and the coexistence region.

The first aspect is the occurrence of instability. In a two-component mixture, one refers to the boundary of this instability as the spinodal, which is a curve. To analyze the local stability of the mixture against phase separation, it is convenient to introduce the so-called Hessian matrix ${\mathit{M}}_{\mathbf{1}}$. The Hessian matrix characterizes the local curvature of the Helmholtz free energy surface. The limit of stability is reached when the matrix has one zero eigenvalue and is otherwise a positive definite (Heidemann and Khalil, 1980). For a mixture with N distinguishable components, the Hessian matrix can be represented by an N × N matrix, referred to as $M_1$ according to

$M_1=\left[\begin{array}{ccc}{2B}_{11}+\frac{1}{c_1}& \cdots & {2B}_{1N}\\ ⋮& \ddots & ⋮\\ {2B}_{1N}& \cdots & {2B}_{NN}+\frac{1}{c_N}\end{array}\right].(5)$

The second aspect to consider is the critical point. In a binary mixture, the critical point is a stable point which lies on the stability limit (spinodal) (Heidemann and Khalil, 1980) and where the phase boundary and spinodal coincide. In mixtures with more than two components, these critical points are also called plait points.

Critical points and plait points are, in general, concentrations at which two or more phases are in equilibrium and become indistinguishable (Heidemann, 1994).

There are two criteria that have to be fulfilled to find the plait points. The first one is det(M₁) = 0, which is the equation for the spinodal. The other criterion is based on the fact that at the critical point, the third derivative of free energy should also be zero. For a multicomponent system, this criterion can be reformulated using Legendre transforms as det(M₂) = 0 (Beegle et al., 1974; Reid and Beegle, 1977), where

$M_2=\left[\begin{array}{ccc}\frac{{\partial \mu }_1}{{\partial n}_1}& \cdots & \frac{{\partial \mu }_n}{{\partial n}_n}\\ ⋮& \ddots & ⋮\\ \frac{{\partial M}_1}{{\partial n}_1}& \cdots & \frac{{\partial M}_1}{{\partial n}_N}\end{array}\right],(6)$

and matrix M₂ is matrix M₁ with one of the rows (in this case the lowest row) being replaced by the partial derivatives of the determinant of matrix M₁. Note that it does not matter which row of the matrix is replaced.

The third aspect is the coexistence region. This is characterized by the equilibrium composition of the phases for N components, which is in turn characterized by a manifold (for two components, this becomes the binodal or coexistence curve). For a system that separates into two different phases, this manifold is found by simultaneously solving Eqs 7, 8:

${\mathrm{\Pi }}^I={\mathrm{\Pi }}^{II},(7)$

${{\mu }_i}^I={{\mu }_i}^{II}\text{with\hspace{0.17em}}i=\mathrm{1,2},\dots ,N.(8)$

For a mixture containing $N$ distinguishable components, which de-mixes into two phases, this yields ($N$ + 1) equations with 2 × $N$ unknowns. In order to solve not only for coexistence concentrations but also for phase volumes, we require extra equations. For the extra set of equations, we build on the fact that no particles are lost and no new particles are created during phase separation and the fact that we assume the total volume does not change. For a system that separates into two separate phases, indicated by I and II, we obtain an extra set of equations (Sturtewagen, 2020):

$\begin{array}{c}{c}_1={\alpha }^I{c}_1^I+\left(1-{\alpha }^I\right)c_1^{II},\\ \cdot \hfill \\ \cdot \hfill \\ \cdot \hfill \\ c_N={\alpha }^I{c}_N^I+\left(1-{\alpha }^I\right)c_N^{II},\end{array}(9)$

where ${\alpha }^I=\frac{V^I}{V^I+V^{II}}$, in which $V^I$ represents the volume of phase $I$.

The solutions for the two-component mixtures are discussed elsewhere (Ersch et al., 2016; Dewi et al., 2020; Dewi et al., 2021; Bot et al., 2021a; Bot et al., 2021b). Eqs. 7–(9) allow for the determination of the concentration of each component in each of the phases along with the phase volume of each phase for any given parent concentration, given that the mixture will separate into two phases.

3 Materials and methods

Size exclusion chromatography with multi-angle laser light scattering (SEC-MALLS) for both polyethylene oxide (PEO35) and dextran (D100) was performed by NIZO food research to obtain the molar mass and size distribution (size distributions are shown in Figure 1). The SEC-MALLS molar mass plots can be found in the Supplementary Materials. The experimental critical point and osmotic virial coefficients for the pure components PEO35 and D100 were taken from Dewi et al. (2020) (see Table I). The cross-virial coefficients used in the calculations are the result of a fitting that includes polydispersity as follows.

FIGURE 1

FIGURE 1. Size distribution of the polydisperse samples that have been used to determine the phase diagram.

Both components exhibited polydispersity (see Figure 1). The size distributions for both components were binned [the number of bins (N = N_PEO + N_D) was N_PEO = 1, 2, and 5 for PEO and N_D = 1, 2, and 10 for dextran]. The obtained radii and fractions were taken as the starting point for the fit of the theoretical polydisperse virial coefficients to the experimentally obtained ones. The cross-interactions between the sub-components of PEO and the sub-components of dextran were also considered the same for all sizes (∆_PD). For each bin size, one can calculate a virial coefficient and the corresponding cross-virial coefficients with the other “components”. Each virial coefficient resulted from a non-additive hard sphere interaction according to Eq. 3 for B_ij, where i can refer to D100 and PEO35 and j can refer to D100 and PEO35. In addition, we assume that for the “components,” the values for ${∆}_{D100}$, ${∆}_{PEO35}$, and ${∆}_{PD}$ are independent of the respective radii.

Subsequently, one calculates a number-averaged virial coefficient for dextran and PEO, B_D100 and B_PEO35, respectively, and the cross-virial coefficient, B_PD, for the size distribution. These calculated expressions still contain the unknown values for the respective parameters, ${∆}_{D100}$, ${∆}_{PEO35}$, and ${∆}_{PD}$. One then determines the values for the non-additivity parameters D_D100, D_PEO35, and D_PD, from fitting the calculated number-averaged expressions to the experimental values of B_D100, B_PEO35, and B_PD. Once the non-additivity parameters are known, we then determine the coexistence relationships and the instability line, and the critical point, and compare these with the ones obtained experimentally by conducting phase diagram studies. One can explore this procedure for the different binning strategies. While using the fitting procedure mentioned above, the sub-component fractions are obtained from the binned-size distribution, while adjusting the sizes of the sub-components within a small range and the non-additivity parameters (∆_PEO35, ∆_D100, and ∆_PD). The values for the virial coefficients and critical points that were used in the fitting are given in Table 1. The obtained matrix of virial coefficients was used to calculate the phase boundary and fractionation of the polydisperse mixture. The resulting phase diagram was compared to the experimental phase diagram from Dewi et al. (2020). A fitting was considered good when the critical point of the polydisperse mixture was on the two-phase boundary and the dilution line through the critical point reached an equal volume at the critical point.

TABLE 1

TABLE 1. Virial coefficients of polyethylene oxide (PEO35) and dextran (D100) and the critical point for their mixture obtained from experiments by Dewi et al. (2020), which were used in the fitting for the polydisperse virial-coefficient matrix.

4 Results and discussion

From Figure 1, we conclude that both PEO and dextran show significant polydispersity. The distribution for PEO is relatively narrow, is slightly bimodal, and has a small tail, while the distribution for dextran has more peaks and a considerably fat tail. Based on the research we did on the effect of polydispersity on the phase behavior of a binary mixture of non-additive hard spheres (Sturtewagen and van der Linden, 2022), we hypothesize that especially the larger components of dextran modify the phase behavior. The osmometric number-averaged virial coefficient of PEO is considerably larger than the number-averaged virial coefficient of dextran (Table 1), even though the radius of the PEO molecules is not larger. This indicates that the repulsive depletion interaction between PEO molecules is larger than the repulsive interaction between dextran molecules.

To determine which effect has the largest influence on the phase behavior (the large size difference between dextran molecules or the higher repulsive interaction between PEO molecules), we have stepwise introduced more polydispersity into our theoretical calculations of the phase diagram. We have compared the calculations to the experimental phase diagram from Dewi et al. (2020). We performed analyses for three cases. First, we introduced polydispersity only for PEO (Supplementary Materials), subsequently we introduced polydispersity only for dextran (Supplementary Materials), and third, we introduced polydispersity for both PEO and dextran (Supplementary Materials). We have compared the different binning strategies as shown in Figure 2; in the main text, we show the result that best fits all the three cases of the experimental coexistence and instability lines. This best fit was obtained using polydispersity of dextran as modeled by the means of 10 bins of equal width, and polyethylene oxide being modeled as monodisperse.

FIGURE 2

FIGURE 2. Coexistence line (black dotted) and instability line (black) for the two-component system of PEO and dextran obtained from solving Eqs 7–9 while substituting the three experimentally determined virial coefficients, as determined by osmometry (see also Dewi et al., 2020). The blue dotted line and the blue line represent the solutions of Eqs.(7)–(9) while using the binning strategy for PEO and dextran, as indicated in the inset of the figure and further described in the text. In fact, PEO is used as monodisperse in the specific strategy depicted here. The circular open symbols represent the experimentally determined one-phase systems and the circular closed symbols represent the two-phase systems.

From Figure 2, we can conclude that the incorporation of polydispersity draws the (apparent) spinodal more toward the horizontal axis and a typical W-form emerges. We note that the term (apparent) spinodal stems from the fact that we plot the concentrations as determined experimentally, effectively hiding the polydispersity of the two polymers, which, if polydispersity would have been included, would, in fact, lead to a multi-dimensional phase diagram instead of a two-dimensional phase diagram.

When subdividing dextran into 10 bins of equal bin width, we capture the details of polydispersity in our fitting. This causes the spinodal to bend toward the vertical axis at higher dextran concentrations, indicating that at higher concentrations of dextran, the mixture can de-mix into two phases of dextran. The binodal of the mixture with this fitting also changes quite drastically. At the dextran-enriched side, the phase boundary first shifts toward higher PEO concentrations and then bends toward the axis with increasing dextran concentrations. The curve has a noticeable nod. This indicates a regime change in demixing. At the PEO-depleted side of the curve, two-phase demixing is with a phase enriched in PEO and a phase enriched in dextran, but it has a phase enriched in the smaller components of dextran and a phase enriched in the larger components of dextran, similar to that described for PEO above. At the PEO-enriched side of the curve, the binodal shifts toward lower concentrations of dextran while asymptotically approaching the horizontal axis. Furthermore, the components with this pairwise interaction can de-mix into more than two phases; however, just as in the previous cases, the concentrations required for multi-phase demixing are unattainable with this particle distribution. The fitting captures more of the experimentally phase-separated samples with low dextran concentrations (filled circles in the plot) compared to the previous fittings. This indicates that our hypothesis that the polydispersity of dextran would be a driving factor in the demixing of these samples is valid.

Döbert et al. (1995) constructed a consistent osmotic virial equation to predict the phase behavior of a different type of polydisperse dextran and monodisperse PEO. They obtained an average molecular weight osmotic virial coefficient for the different polymer chains and compared their results to the predictions, assuming a binary monodisperse system. They report a better fit when polydispersity is taken into account. When predicting phase separation for mixtures close to the phase boundary, they also report that the binodal shifts closer to the axis in the dextran-depleted side of the binodal. They report a strong influence of the volume fraction of the phase enriched in dextran on the shape of the phase boundary and the fractionation of polydisperse dextran.

Regarding the validity of the virial approach up to the second order, we like to remark that we determine the second virial coefficients from the slope of the reduced osmotic pressure, $\left(\mathrm{\Pi }/RTc\right)$ as a function of $c$, in the limit of $c$ approaching zero. The contributions of the higher order virial coefficients can be ignored for our system since the slope was found to be constant in the concentration range that was used. We note that the osmotic pressures were determined in the concentration range of 0.02–0.2 mol/m³, which is well below the critical concentrations: (PEO35; D100) = (0.4; 1.06) (Dewi et al., 2020). Interpreting the osmotic pressure experiments in terms of the polymers being modeled by non-additive hard spheres is justified as follows. The (estimated) overlap concentration of the polymers, c^*, for this system, lies in the order of 0.5 mol/m³ when assuming a radius of gyration of 10 nm. This radius of gyration is taken from the size distribution in Figure 1. The measurements for the virial coefficients are conducted well below this overlap concentration. The critical concentration is in the same range as the overlap concentration, but it is noted in support of the second-order approximation that the results for the critical point, as deduced from substituting the experimental values for the virial coefficients in the theoretical second-order model, yield a consistent picture with the experimental critical point. We, therefore, expect that the second-order virial approach gives reasonable predictions and refers to a more elaborate discussion on this topic in Sturtewagen and van der Linden (2023).

Predicting the phase behavior of multiple components has been of interest in the literature for some time. One approach that has been put forward recently builds on a convenient parameterization that leads to a set of equations that are much more easily analyzed than applying numerical methods to solve all equations (Bot et al., 2023). We briefly summarize this work in the following paragraph.

First, consider for simplicity the binodal of a two-component system. Two coexisting phases, named I and II, are represented by two molar concentration coordinates $\left(c_1^I\text{\hspace{0.17em}and\hspace{0.17em}}{c}_2^I\right)$ and $\left(c_1^{II}\text{\hspace{0.17em}and\hspace{0.17em}}{c}_2^{II}\right)$ on the binodal, which are connected by a so-called tie-line. Introducing the parameter $S_{m,21}=-\left(c_2^{II}-c_2^I\right)/\left(c_1^{II}-c_1^I\right)$, corresponding to (minus) the tie-line slope, allows for rewriting the coexistence equations. Introducing the parameter $S_{m,21}$ may seem like a step back, since it adds a fourth equation to the coexistence equations. However, two of these equations have analytical solutions in terms of the Lambert W function (Corless et al., 1996), and as a result, the original four-variable problem is reduced to a more succinct problem in the two concentration coordinates and the parameter $S_{m,21}$. From a physical perspective, the choice of (minus) the tie-line slope as a parameter seems to be a defendable choice, as the tie-line contains important information on the physics of the problem. For real arguments, the Lambert W function has zero, one, or two solutions, corresponding to, respectively, the isotropic mixing, location of the critical points, and binary phase separation. In case of two solutions, the solution for each component in phase I and II is located either on the W₋₁-branch or W₀ branch of the Lambert W function. Surprisingly, this approach can be generalized to many components, i.e., for higher values than N = 2. There are N² parameters $S_{m,ij}=-\left(c_j^{II}-c_j^I\right)/\left(c_i^{II}-c_i^I\right)$, of which (N − 1) are free variables and $S_{m,ii}=-1$ is fixed (for components i and j). This reduces the number of variables in the problem from 2N concentration variables to N concentration variables plus (N − 1) variables $S_{m,ij}$. For more details, we refer to a recent publication of Bot et al. (2023). The results in Table 2 present a scheme to calculate the phase diagram using the abovementioned parametrization. Previously, results have been published for N = 2 (Bot et al., 2021a; Bot et al., 2021b) and N = 3 (Bot and Venema, 2023). A comparison between the calculations in terms of the parameters $S_{m,ij}$ (Table 2) and the numerical results for polydisperse two-component mixtures of PEO and dextran, which are presented in Section 4.2, remains to be evaluated. The calculation of the binodal manifolds in particular pose numerical challenges because of the large number of coupled non-linear algebraic equations. The calculation of the critical and spinodal manifolds presents fewer complications. For an arbitrary, large number of components N, it can be considered complementary to another interesting method, referred to in the literature as the random matrix theory (RMT) (Sear and Cuesta, 2003; Thewes et al., 2023).

TABLE 2

TABLE 2. Expressions for the spinodal, critical, and binodal manifolds for an arbitrary number of components N ≥ 2 in a common solvent, in terms of the parameters ${\mathrm{S}}_{\mathrm{m},\text{ij}}$ and the Lambert W function, where i = 1, …, N, and z can be represented as “sp”, “c”, or “m” (adapted from Bot et al., 2023).

This random matrix theory has been put forward in the context of phase behavior two decades ago by Sear and Cuesta (2003) as an alternative approach to handle the complexity of mixtures with many components, as a utilization of the general work by Wigner (1967). Setting details aside, this method allows for calculating the spinodal curve of complex mixtures using only averages and standard deviations of the second virial coefficients of the molecules in the mixture. This simplifies the calculations considerably. The RMT uses essentially the same expression for free energy as shown in Eq. 1. Additionally, it makes use of the large number of components and accordingly the large number of second virial coefficients. Instead of finding solutions for the spinodal manifold (with the result shown in the left column of Table 2), Sear and Cuesta (2003) calculated a point on the spinodal along the line $c_1=c_2=\dots =c_N$, where $c_i$ is the molar concentration of component i, for a system of many components N. Note that because the number of components N is large, the individual concentrations of the components are very small. This approach was recently extended to the arbitrary concentrations of the components by Thewes et al. (2023). The second virial coefficients $B_{ij}$ are assumed to form a (symmetric) matrix, composed of diagonal elements that can be randomly chosen due to the large number of components present, allowing randomization in this choice, with an average value b’ and a standard deviation σ and off-diagonal elements with an average value b and also a standard deviation σ. The virial coefficients satisfy the criterion of statistical independence, and their absolute value is smaller than a finite value. These assumptions ensure that the virial coefficients, i.e., the matrix components, can be randomly chosen according to a predefined distribution. This considerably simplifies many aspects of the calculations of the spinodal manifold of such mixtures. The premise of the RMT is that if the random matrix is large enough, the actual values of $B_{ij}$ typically do not matter anymore and only their average and standard deviations are of importance, and for some aspects, even the details of the distribution itself do not play a role anymore. It is not completely clear, however, which simplifications in the calculations are driven by mathematical convenience and which ones are driven by physical relevance, thus yielding some words of caution while applying the RMT approach, as was also mentioned by Jacobs (2023). The exact analytical expressions given in Table 2 allow us to separate the mathematical aspects of the RMT approach from the physical aspects in the calculation. In particular, the results in Table 2 encompass the critical point, the coexistence manifold, and the spinodal manifold, whereas the RMT only yields the spinodal manifold, often with more severe assumptions regarding the concentrations of the different components being all equal. For the present article, a direct comparison between the predictions of the RMT and Table 2 is out of scope.

5 Conclusion

The equilibrium phase behavior of multicomponent systems is relevant for understanding gelation of such mixtures. Controlling the gel properties can be attained via controlling the structures being formed during the non-equilibrium processes taking place during gel formation, while crossing equilibrium phase boundaries.

We have described the theoretical (numerical) work on predicting the equilibrium phase behavior of a multicomponent system consisting of polydisperse polymers and compared this to the experimental data on the phase diagram. Taking into account the polydispersity helps to more accurately predict the particular form of the spinodal and binodal.

Polydispersity plays an important role in the phase behavior of the polydisperse polymers PEO and dextran. The components with a larger-than-average molecular weight govern the transition between one and two phases close to the phase boundary in their respective depleted concentration ranges. This causes a decrease in miscibility and a shift of the phase boundary to lower concentrations. This causes drastic changes to the shape of the phase boundary. When both components are polydisperse, the phase boundary drastically changes shape, and changes from a U-shape to a W-shape. It is not only the phase boundary that changes shape, but the spinodal curve also has different boundaries depending on the polydispersity. Even though multi-phase separations for mixtures with the fitted interactions are possible, the existing particle size distributions of PEO and dextran make concentrations resulting in multi-phase systems unattainable.

Upon demixing, the distribution of polydisperse components changes in each phase (Supplementary Materials). This fractionation is dependent on the parent distribution, the pair-wise interaction between the components of the same type, the pair-wise interaction of the components of a different type, and the concentration of both components in the parent mixture.

Our method of incorporating polydispersity allows for a more precise prediction of the phase boundary compared to assuming monodispersity, especially in the metastable region. Next to that, our method allows for prediction of the concentration and fractionation of each component in each phase depending on the parent concentration and the volume fraction of the said phases.

The best fit with the available data to the experimental data was when dextran was polydisperse and PEO was monodisperse.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

LS: conceptualization, investigation, methodology, manuscript writing—review and editing, data curation, software, and visualization. BPCD: investigation, methodology, visualization, manuscript writing—review and editing, and funding acquisition. AB: conceptualization, formal analysis, investigation, methodology, and manuscript writing—review and editing. PV: conceptualization, formal analysis, investigation, methodology, software, supervision, and manuscript writing—review and editing. EvL: conceptualization, funding acquisition, investigation, methodology, project administration, supervision, manuscript writing—original draft, review, and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. BPCD acknowledges funding from the Indonesia Endowment Fund for Education (LPDP—Lembaga Pengelola Dana Pendidikan) scholarship, Ministry of Finance, The Republic of Indonesia.

Conflict of interest

Author AB was employed by Unilever.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, editors, and reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsfm.2023.1328180/full#supplementary-material

References

Albertsson, P.-Å. (1970). Partition of cell particles and macromolecules in polymer two-phase systems. Adv. Protein Chem. 24, 309–341. doi:10.1016/s0065-3233(08)60244-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bagnani, M., Nyström, G., De Michele, C., and Mezzenga, R. (2019). Amyloid fibrils length controls shape and structure of nematic and cholesteric tactoids. ACS Nano 13 (1), 591–600. doi:10.1021/acsnano.8b07557

PubMed Abstract | CrossRef Full Text | Google Scholar

Beegle, B. L., Modell, M., and Reid, R. C. (1974). Thermodynamic stability criterion for pure substances and mixtures. AIChE J. 20 (6), 1200–1206. doi:10.1002/aic.690200621

CrossRef Full Text | Google Scholar

Bot, A., Dewi, B. P. C., and Venema, P. (2021a). Addition to “phase-separating binary polymer mixtures: the degeneracy of the virial coefficients and their extraction from phase diagrams”. ACS Omega 6 (30), 20086–20087. doi:10.1021/acsomega.1c03331

PubMed Abstract | CrossRef Full Text | Google Scholar

Bot, A., Dewi, B. P. C., and Venema, P. (2021b). Phase-separating binary polymer mixtures: the degeneracy of the virial coefficients and their extraction from phase diagrams. ACS Omega 6 (11), 7862–7878. doi:10.1021/acsomega.1c00450

PubMed Abstract | CrossRef Full Text | Google Scholar

Bot, A., van der Linden, E., and Venema, P. (2023). Phase separation in complex mixtures with many components: analytical expressions for spinodal manifolds. Arxiv 2306.17600. Available at: https://doi.org/10.48550/arXiv.2306.17600.

Google Scholar

Bot, A., and Venema, P. (2023). Phase behavior of ternary polymer mixtures in a common solvent. ACS Omega 8 (31), 28387–28408. doi:10.1021/acsomega.3c02604

PubMed Abstract | CrossRef Full Text | Google Scholar

Chatellier, J. Y., Durand, D., and Emery, J. R. (1985). Critical helix content in gelatin gels. Int. J. Biol. Macromol. 7 (5), 311–314. doi:10.1016/0141-8130(85)90030-3

CrossRef Full Text | Google Scholar

Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J., and Knuth, D. E. (1996). On the LambertW function. Adv. Comput. Math. 5 (1), 329–359. doi:10.1007/bf02124750

CrossRef Full Text | Google Scholar

Dewi, B. P. C., van der Linden, E., Bot, A., and Venema, P. (2020). Second order virial coefficients from phase diagrams. Food Hydrocoll. 101, 105546. doi:10.1016/j.foodhyd.2019.105546

CrossRef Full Text | Google Scholar

Dewi, B. P. C., van der Linden, E., Bot, A., and Venema, P. (2021). Corrigendum to “second order virial coefficients from phase diagrams.”[Food Hydrocolloids 101 (2020) 105546]. Food Hydrocoll. 112, 106324. doi:10.1016/j.foodhyd.2020.106324

CrossRef Full Text | Google Scholar

Dijkstra, M., Brader, J. M., and Evans, R. (1999). Phase behaviour and structure of model colloid-polymer mixtures. J. Phys. Condens. Matter 11, 10079–10106. doi:10.1088/0953-8984/11/50/304

CrossRef Full Text | Google Scholar

Döbert, F., Pfennig, A., and Stumpf, M. (1995). Derivation of the consistent osmotic virial equation and its application to aqueous poly(ethylene glycol)–dextran two-phase systems. Macromolecules 28 (23), 7860–7868. doi:10.1021/ma00127a037

CrossRef Full Text | Google Scholar

Dos Santos, T. E., Loh, W., and Pessôa-Filho, P. D. A. (2004). Phase equilibrium in aqueous two-phase systems containing ethylene oxide-propylene oxide block copolymers and dextran. Fluid Phase Equilibria 218 (2), 221–228. doi:10.1016/j.fluid.2004.01.001

CrossRef Full Text | Google Scholar

Edelman, M. W., Tromp, R. H., and van der Linden, E. (2003). Phase-separation-induced fractionation in molar mass in aqueous mixtures of gelatin and dextran. Phys. Rev. E 67 (2), 021404. doi:10.1103/physreve.67.021404

PubMed Abstract | CrossRef Full Text | Google Scholar

Ersch, C., Meijvogel, L. L. C., van der Linden, E., Martin, A., and Venema, P. (2016). Interactions in protein mixtures. Part I: second virial coefficients from osmometry. Food Hydrocoll. 52, 982–990. doi:10.1016/j.foodhyd.2015.07.020

CrossRef Full Text | Google Scholar

Gaube, J., Höchemer, R., Keil, B., and Pfennig, A. (1993). Polydispersity effects in the system poly(ethylene glycol) + dextran + water. J. Chem. Eng. Data 38 (2), 207–210. doi:10.1021/je00010a005

CrossRef Full Text | Google Scholar

Heidemann, R. A. (1994). The classical theory of critical points. Dordrecht: Springer, 39–64.

CrossRef Full Text | Google Scholar

Heidemann, R. A., and Khalil, A. M. (1980). The calculation of critical points. AIChE J. 26 (5), 769–779. doi:10.1002/aic.690260510

CrossRef Full Text | Google Scholar

Jacobs, W. M. (2023). Theory and simulation of multiphase coexistence in biomolecular mixtures. J. Chem. Theory Comput. 19, 3429–3445. doi:10.1021/acs.jctc.3c00198

PubMed Abstract | CrossRef Full Text | Google Scholar

Johansson, G., and Walter, H. (1999). Partitioning and concentrating biomaterials in aqueous phase systems. Int. Rev. Cytol. 192 (1896), 33–60. doi:10.1016/s0074-7696(08)60521-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Johansson, H. O., Karlström, G., Tjerneld, F., and Haynes, C. A. (1998). Driving forces for phase separation and partitioning in aqueous two-phase systems. J. Chromatogr. B Biomed. Appl. 711 (1-2), 3–17. doi:10.1016/s0378-4347(97)00585-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Joly-Duhamel, C., Hellio, D., Ajdari, A., and Djabourov, M. (2002a). All gelatin networks: 2. The master curve for elasticity. Langmuir. 18 (19), 7158–7166. doi:10.1021/la020190m

CrossRef Full Text | Google Scholar

Joly-Duhamel, C., Hellio, D., and Djabourov, M. (2002b). All gelatin networks: 1. Biodiversity and physical chemistry. Langmuir. 18 (19), 7208–7217. doi:10.1021/la020189n

CrossRef Full Text | Google Scholar

Kang, C. H., and Sandler, S. I. (1988). Effects of polydispersivity on the phase behavior of the aqueous two-phase polymer systems. Macromolecules 21 (10), 3088–3095. doi:10.1021/ma00188a029

CrossRef Full Text | Google Scholar

Larsson, M., and Mattiasson, B. (1988). Characterization of aqueous two-phase systems based on polydisperse phase forming polymers: enzymatic hydrolysis of starch in a PEG–starch aqueous two-phase system. Biotechnol. Bioeng. 31 (9), 979–983. doi:10.1002/bit.260310910

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsson, M., and Mattiasson, B. (1992). Molecular weight distribution analysis of water-soluble polymers using aqueous two-phase systems. Ann. N. Y. Acad. Sci. 672 (1), 649–656. doi:10.1111/j.1749-6632.1992.tb32741.x

CrossRef Full Text | Google Scholar

McMillan, W. G., and Mayer, J. E. (1945). The statistical thermodynamics of multicomponent systems. J. Chem. Phys. 13 (7), 276–305. doi:10.1063/1.1724036

CrossRef Full Text | Google Scholar

Nguyen, M., and Vaikuntanathan, S. (2016). Design principles for nonequilibrium self-assembly. Proc. Natl. Acad. Sci. U. S. A. 113 (50), 14231–14236. doi:10.1073/pnas.1609983113

PubMed Abstract | CrossRef Full Text | Google Scholar

Normand, V., Muller, S., Ravey, J.-C., and Parker, A. (2000). Gelation kinetics of gelatin: A master curve and network modeling. Macromolecules. 33 (3), 1063–1071. doi:10.1021/ma9909455

CrossRef Full Text | Google Scholar

Normand, V., and Parker, A. (2003). “Scaling in the dynamics of gelatin gels,” in Proceedings of the third International Symposium on Food Rheology and Structure Lappersdorf, Germany, 2003 (Kerschensteiner Verlag).

Google Scholar

Poulin, P., Bibette, J., and Weitz, D. A. (1999). From colloidal aggregation to spinodal decomposition in sticky emulsions. Eur. Phys. J. B. Condensed. Matter. Complex. Systems. 7 (2), 277–281. doi:10.1007/s100510050614

CrossRef Full Text | Google Scholar

Reid, R. C., and Beegle, B. L. (1977). Critical point criteria in legendre transform notation. AIChE J. 23 (5), 726–732. doi:10.1002/aic.690230515

CrossRef Full Text | Google Scholar

Sear, R. P., and Cuesta, J. A. (2003). Instabilities in complex mixtures with a large number of components. Phys. Rev. Lett. 91 (24), 245701. doi:10.1103/physrevlett.91.245701

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturtewagen, L. (2020). Predicting phase behavior of multi-component and polydisperse aqueous mixtures using a virial approach. PhD Thesis. Wageningen: Wageningen University.

Google Scholar

Sturtewagen, L., and van der Linden, E. (2021). Effects of polydispersity on the phase behavior of additive hard spheres in solution. Molecules 26 (6), 1543. doi:10.3390/molecules26061543

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturtewagen, L., and van der Linden, E. (2022). Towards predicting partitioning of Enzymes between macromolecular phases: effects of polydispersity on the phase behavior of nonadditive hard spheres in solution. Molecules 27 (19), 6354. doi:10.3390/molecules27196354

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturtewagen, L., and van der Linden, E. (2023). Ternary mixtures of hard spheres and their multiple separated phases. Molecules 28 (23), 7817. doi:10.3390/molecules28237817

PubMed Abstract | CrossRef Full Text | Google Scholar

Thewes, F. C., Krüger, M., and Sollich, P. (2023). Composition dependent instabilities in mixtures with many components. Phys. Rev. Lett. 131 (5), 058401. doi:10.1103/physrevlett.131.058401

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Linden, E. (2012). From peptides and proteins to micro-structure mechanics and rheological properties of fibril systems. Food Hydrocoll. 26 (2), 421–426. doi:10.1016/j.foodhyd.2010.11.019

CrossRef Full Text | Google Scholar

van der Linden, E., and Parker, A. (2005). Elasticity due to semiflexible protein assemblies near the critical gel concentration and beyond. Langmuir. 21 (21), 9792–9794. doi:10.1021/la051312o

PubMed Abstract | CrossRef Full Text | Google Scholar

Wigner, E. P. (1967). Random matrices in physics. SIAM Rev. 9 (1), 1–23. doi:10.1137/1009001

CrossRef Full Text | Google Scholar

Zhao, Z., Li, Q., Ji, X., Dimova, R., Lipowsky, R., and Liu, Y. (2016). Molar mass fractionation in aqueous two-phase polymer solutions of dextran and poly(ethylene glycol). J. Chromatogr. A 1452, 107–115. doi:10.1016/j.chroma.2016.04.075

PubMed Abstract | CrossRef Full Text | Google Scholar