Phase transitions in biological systems with many components

TL;DR

This paper demonstrates through numerical evidence that biological mixtures with many components can spontaneously phase separate into distinct regions, with the process being robust and tunable, facilitating spatial organization in cells.

Contribution

It reveals that multicomponent biological mixtures naturally tend to undergo phase separation due to general phase diagram features, highlighting robustness and tunability of cellular compartmentalization.

Findings

Demixed domains form when interaction variance exceeds a threshold.

Multiple phases can coexist under similar conditions, tunable by composition.

Phase separation remains robust with increasing component number.

Abstract

Biological mixtures such as the cytosol may consist of thousands of distinct components. There is now a substantial body of evidence showing that, under physiological conditions, intracellular mixtures can phase separate into spatially distinct regions with differing compositions. In this paper we present numerical evidence indicating that such spontaneous compartmentalization exploits general features of the phase diagram of a multicomponent biomolecular mixture. In particular, we show that demixed domains are likely to segregate when the variance in the inter-molecular interaction strengths exceeds a well-defined threshold. Multiple distinct phases are likely to become stable under very similar conditions, which can then be tuned to achieve multiphase coexistence. As a result, only minor adjustments to the composition of the cytosol or the strengths of the inter-molecular interactions are needed to regulate the formation of different domains with specific compositions, implying that phase separation is a robust mechanism for creating spatial organization. We further predict that this functionality is only weakly affected by increasing the number of components in the system. Our model therefore suggests that, for purely physico-chemical reasons, biological mixtures are naturally poised to undergo a small number of demixing phase transitions.

Equations

$$\theta \equiv {\cos }^{-1}\left[\frac{{\mathbf{\rho }}^{(0)}\cdot \left({\mathbf{\rho }}^{(2)}-{\mathbf{\rho }}^{(1)}\right)}{\Vert {\mathbf{\rho }}^{(0)}\Vert \Vert {\mathbf{\rho }}^{(2)}-{\mathbf{\rho }}^{(1)}\Vert }\right].$$

$$U=-\sum _{(u,v)\in \{\text{n.n.}\}}\left[\sum _{i=1}^N\sum _{j=1}^N{\delta }_{C_{u}i}{ϵ}_{ij}{\delta }_{C_{v}j}\right].$$

$$\beta F=\sum _i{\rho }_i\ln {\rho }_i+(1-\varphi )\ln (1-\varphi )-\frac{\beta z}{2}\sum _{ij}{\rho }_i{ϵ}_{ij}{\rho }_j,$$

$$\frac{\beta \Delta F}{\varphi }=\left(\sum _i{x}_i\ln x_i+\ln N\right)-\frac{\beta z\varphi }{2}\sum _{ij}{x}_i\Delta {ϵ}_{ij}{x}_j,$$

$$\frac{\Delta F(n)}{{\beta }^{-1}\varphi }=\ln \left(\frac{N}{n}\right)-{\beta }^{\prime }\max {\left[\frac{1}{n^2}\sum _{ij}{\tau }_i^{(n)}\Delta {ϵ}_{ij}{\upsilon }_j^{(n)}\right]}_{\{{\mathbf{\tau }}^{(\mathbf{n})},{\mathbf{\upsilon }}^{(\mathbf{n})}\}},$$

Full text

Introduction

Biological systems carry out complex chemical reactions involving a large number of interacting components. However, living cells should not be viewed as well-mixed reaction vessels. Rather, many of the molecular species in a cell exist as phase-separated domains (1; 2), and there is considerable evidence that phase separation plays a functional role. Examples include the formation of cytoplasmic granules (3; 4; 5; 6), nucleoli (7; 8) and amorphous clusters of proteins involved in signaling pathways (9). These structures typically appear as liquid-like droplets that are selectively enriched in some components, resulting in spontaneous compartmentalization without the use of membranes (10). In some cases, more than two compositionally distinct domains, such as the concentric droplets observed in the nucleolus (11), can form simultaneously. Similarly, in some lipid membranes, interactions between embedded components are believed to drive the formation of phase-separated rafts (12; 13; 14). For many of these examples, there is strong evidence that phase separation is essentially an equilibrium phenomenon: active processes, such as ATP hydrolysis, are not necessary to observe spatial segregation (15; 16).

If the formation of compositionally distinct domains is important for the normal functioning of a cell, it is logical to ask what conditions are required to achieve a heterogeneous spatial distribution of components. To address this question, we first propose three general criteria for ‘biologically functional’ phase separation. For spontaneous compartmentalization to be useful, each domain should be enriched in a relatively small number of functionally related components. In order to regulate the formation of domains, stabilizing a phase-separated region should only require small adjustments to the abundances and inter-molecular interactions that contribute to the segregated phase. Finally, the desired phase behavior should be robust to fluctuations in the interactions between functionally unrelated components.

To understand the physical requirements for meeting these criteria, we consider what happens in a multicomponent mixture with random pairwise interactions between $N$ distinct components. Adopting a model first proposed by Sear and Cuesta (17), we treat the interactions between the molecules in the mixture as pairwise additive and isotropic, i.e., independent of the relative orientations of the molecules (18). As a result, a single random variable can be used to describe the average strength of the interaction between any pair of components. Assuming no detailed knowledge of the molecular interactions, we draw each of the ${N(N+1)/2}$ random interactions independently from the same distribution (that we shall specify later).

We find that this random-mixture model generates phase diagrams that meet our criteria for biological functionality under rather general conditions. The phase diagrams of random mixtures fall into two classes, depending on the distribution of interactions and the number of components: one supporting demixed phases that are enriched in a relatively small number of components, and a much simpler case where two phases with nearly identical compositions coexist. As predicted by a previous simulation study (19), we demonstrate that increasing the number of components tends to favor the simpler binary phase behavior. However, we find that demixed phases appear suddenly at a threshold that is primarily determined by the variance of the inter-molecular interactions and scales weakly with the number of components. Under the conditions where demixing occurs, we find that a limited number of phases, each enriched in a different subset of components, tend to appear at very similar chemical potentials. This observation implies that, with small changes to the component concentrations, it is possible to turn individual phases on and off and to stabilize multiple phases simultaneously; at the same time, random concentration fluctuations are unlikely to drive demixing into phases with arbitrary compositions. Our model therefore suggests that biologically functional phase separation is easily achieved in multicomponent mixtures and that the selection of specific phases can be accomplished by locally tuning the inter-molecular interactions. We provide an intuitive explanation for these results and then discuss the implications for biological mixtures.

Theory and Monte Carlo Simulations

Phase diagrams in many dimensions

Classical thermodynamics predicts that in a mixture with $N$ components, it is possible to observe as many as ${N+2}$ coexisting equilibrium phases (20). The precise number of coexisting phases that appear in any particular system depends on the interactions between the components. In what follows, we shall use an ‘implicit solvent’ picture in which the solvent is not considered to be a component but is merely the medium in which the various molecular species move. Moreover, we shall assume that both the temperature and pressure are constant, although, in general, the compositions of the coexisting phases are both temperature and pressure-dependent.

To illustrate some of the features of multicomponent phase separation, we sketch three possible phase diagrams for a two-component mixture, i.e., two solutes plus solvent, in Figure 1a–c. Because the components may have dissimilar pairwise interactions, phase separation can drive the formation of phases with differing compositions. It is convenient to characterize an $N$-component mixture by the concentration vector ${\bm{\rho}\equiv\{\rho_{1},\rho_{2},\cdots\!\,,\rho_{N}\}}$, where the total concentration of a phase is ${\phi\equiv\sum_{i=1}^{N}\rho_{i}}$. We use dimensionless units for the component concentrations so that all concentration vectors lie within a unit (${N+1}$)-dimensional simplex, such that ${0\leq\phi\leq 1}$. The composition of a phase refers to the mole fraction of each component, ${x_{i}\equiv\rho_{i}/\phi}$. Consequently, the composition describes the direction (but not the magnitude) of $\bm{\rho}$.

During a phase transition, both the magnitude and the orientation of $\bm{\rho}$ can change. Let us imagine that, prior to phase separation, the initially homogeneous mixture is described by the parent concentration vector $\bm{\rho^{(0)}\!}$. After phase separation, coexistence is established between phases with concentration vectors $\bm{\rho^{(1)}}$ and $\bm{\rho^{(2)}}$ and associated total concentrations $\phi_{1}<\phi_{2}$. In the phase diagrams shown in Figure 1a–c, the shaded areas denote coexistence regions, and the circles indicate pairs of daughter phases $\bm{\rho^{(1)}}$ and $\bm{\rho^{(2)}}$; the parent concentration vector lies on the dashed tie line connecting the two daughter phases. To characterize the nature of a particular phase transition, it is useful to define an ‘angle of phase separation,’ the angle between the parent concentration vector and the tie line connecting the two coexisting phases:

[TABLE]

This angle serves as an order parameter that indicates the similarity in the compositions of the parent phase and the coexisting daughter phases.

Multicomponent phase separation can result in phases with equal compositions (Figure 1a), the formation of a high-concentration phase that is enriched in one component (Figure 1b), or some combination thereof. In the case where ${\theta\rightarrow 0}$, the coexisting phases differ only in their total concentrations. This scenario is referred to as a condensation phase transition. Alternatively, if phase separation at an initially low concentration is driven by the demixing of a single component, then ${\theta\rightarrow\theta_{N}\equiv\cos^{-1}(N^{-1/2}})$. In principle, in mixtures with more than two components, intermediate phase behavior can occur if $\theta$ takes a value between 0 and $\theta_{N}$; this scenario implies the selective phase separation of many, but not all, components. Figure 1c also shows a three-phase coexistence region, where the mixture phase separates into three coexisting phases with concentrations corresponding to the vertices of the dark shaded triangle. For parent concentrations lying within this region, multiple demixing transitions occur at the same temperature and pressure. However, only one phase transition is typically encountered at the boundary of the homogeneous phase, unless the interaction matrix is degenerate.

Characterization of phase diagrams via simulation

Simulations allow us to study the nature of phase separation in the limit of very many components. Unlike the canonical phase diagrams shown in Figure 1a–c, we perform simulations in the grand-canonical ensemble, where the chemical potential of each component is held fixed and the number of particles of each type is allowed to fluctuate. The grand-canonical ensemble is ideal for studying phase transitions because it allows pure phases to occupy the entire simulation volume. This means that we can observe coexistence without forming thermodynamically unfavorable interfaces between the pure phases.

Because the component concentrations fluctuate, a grand-canonical simulation samples a free-energy landscape in which low-free-energy basins correspond to the phases in a canonical phase diagram. For each representative phase diagram shown in Figure 1a–c, a corresponding free-energy landscape is shown below in Figure 1d–f. In each case, the illustration shows the landscape for the component chemical potentials that correspond to a parent concentration vector on the highlighted tie line. Importantly for what follows, simulations in the grand-canonical ensemble reveal metastable free-energy basins even when the corresponding phase is unstable in the canonical phase diagram. For example, the metastable basin in the asymmetric free-energy landscape shown in Figure 1f corresponds to a third phase that must be stable at a nearby set of chemical potentials, indicating that the tie line on the canonical phase diagram in Figure 1c is close to a three-phase coexistence region.

In order to to determine whether a matrix of inter-molecular interactions supports multiphase coexistence, we follow a simulation strategy that allows us to characterize the canonical phase diagram with a limited number of grand-canonical simulations. We first assume that all components are present in equal amounts in the parent phase. Then, starting from the low-concentration homogeneous phase, we locate the lowest-concentration phase boundary that intersects the parent composition vector, i.e., ${\bm{\rho^{(1)}}=\bm{\rho^{(0)}}\!}$. This concentration vector is commonly called the cloud point for a given parent composition and is indicated by the filled circles in Figure 1a–c. Lastly, we sample the free-energy landscape and identify metastable free-energy basins at the cloud-point chemical potentials to determine the proximity of this phase boundary to a multiphase coexistence region. These two sets of simulations are described below.

Monte Carlo simulations of binary phase coexistence

To locate the cloud point associated with an equimolar parent composition, we carried out grand-canonical Monte Carlo (GCMC) simulations (21) of a three-dimensional lattice model. In this model, lattice sites may either be vacant, representing the non-interacting solvent, or occupied by one of the $N$ components. The energy, $U$, of a lattice configuration is calculated by summing the interactions between all pairs of particles on nearest-neighbor (n.n.) lattice sites:

[TABLE]

The tuples ${(u,v)}$ run over all nearest-neighbor pairs, ${C_{u}}$ is the component index of the particle at lattice site $u$ (${C_{u}=0}$ if the site is vacant), and $\delta$ is the Kronecker delta.

We employed multicanonical biasing (22) to facilitate rapid crossing of the free-energy barriers that separate low and high-concentration lattice configurations. In order to locate the cloud point in these simulations, we tuned the component chemical potentials to achieve equal free energies in the low-concentration phase and the most stable high-concentration phase, ${F^{(1)}=F^{(2)}\!}$, while satisfying the imposed composition constraint on the low-concentration phase, ${\bm{x^{(1)}}=\bm{x^{(0)}}}$. This strategy has been shown to minimize finite-size effects (23). We note that this approach does not constrain the simulation to sample only from the two phases of interest if other free-energy basins are present. For computational efficiency, we used a ${L\times L\times L}$ cubic lattice with periodic boundaries and ${L=6}$. Complete details of the simulation method are provided in Ref. 19.

Free-energy landscapes near multiphase coexistence

For each realization of the interaction matrix, we generated approximately 1000 random lattice configurations, with concentration vectors chosen uniformly from the unit simplex of component concentrations. We then allowed each initial configuration to evolve via the unbiased GCMC algorithm, with the component chemical potentials fixed at the previously determined cloud point, until the concentration fluctuations stabilized. These GCMC trajectories tend to travel ‘downhill’ on the free-energy landscape, and, without the multicanonical biasing that was used in the previous set of simulations, the probability of escaping from a free-energy basin is extremely small. The endpoints of these trajectories therefore lie close to the minima of the free-energy basins, which correspond to stable or metastable thermodynamic phases. We clustered these endpoints in order to determine the number of distinct free-energy basins, and then we found the mean concentration vector of each of these (meta)stable phases by averaging the endpoint concentration vectors within each basin.

Results

Simulations of condensation and demixing in mixtures with random pairwise interactions

To study an ensemble of mixtures with random interactions, we generated many independent realizations of the pairwise interaction matrices ${\bm{\epsilon}\equiv\{\epsilon_{ij}\}}$. The interactions between pairs of components were assumed to be independent random variables drawn from a Gaussian distribution with mean $\bar{\epsilon}$ and variance $\sigma^{2}$. (We shall discuss the effects of correlated interactions in the context of the mean-field theory presented below.) We limited the number of components to 64 for computational tractability and then chose the variance of the Gaussian distributions in order to observe both the demixing and condensation limits of the phase behavior. The mean interaction strength was fixed at ${\bar{\epsilon}=1.07\epsilon_{\text{c}}}$, where ${-\epsilon_{\text{c}}=-0.87k_{\text{B}}T}$ is the critical bond energy of the one-component lattice gas, $k_{\text{B}}$ is the Boltzmann constant and $T$ is the absolute temperature. We shall discuss the relevance of the critical temperature and the rationale for choosing this mean interaction strength in the mean-field section below.

For each realization of ${\bm{\epsilon}}$, we used our simulation strategy both to calculate the component chemical potentials at the lowest-concentration phase boundary and then to determine the number of free-energy basins at this point. We also computed the average composition of the coexisting high-concentration phase at this set of chemical potentials in order to calculate $\theta$. As noted above, it is highly unusual to observe multiple coincident transitions precisely at this first phase boundary. However, when we do observe more than two free-energy basins, multiphase coexistence must occur at very similar component chemical potentials. In such a case, a stronger driving force — either higher concentrations or stronger average interactions — is necessary to push the mixture into the multiphase region. This scenario is consistent with our expectations for functional multiphase coexistence, because phase transitions that are close in chemical-potential space may be easily manipulated by small changes to the component abundances and inter-molecular interactions.

Bimodal distribution of phase behaviors

Our simulations show that random interactions result in two distinct types of phase behavior: either a handful of components will demix, or the mixture will condense into two phases with nearly identical compositions. The resulting bimodal distribution of the angle of phase separation can be seen in Figure 2a–b, where we have included all high-concentration phases that were identified from the free-energy landscapes of dozens of random mixtures with 32 or 64 components. Low values of $\theta$ correspond to condensation phase transitions. In this limit, a multicomponent mixture behaves as if all the components were alike. In contrast, the higher-$\theta$ peak corresponds to demixing phase transitions; in this case, the maximum angle of phase separation is bounded by $\theta_{N}$. Importantly, in simulations with 8 or more components, we rarely find phase behavior that is intermediate between condensation and demixing.

The average phase behavior of a random mixture ensemble depends on both the width of the distribution of inter-molecular interactions, $\sigma$, and the number of components, $N$. Condensation transitions dominate for large $N$ and small $\sigma$, while demixing transitions are more likely in mixtures with fewer components and a broader distribution of interactions. In fact, the average angle of phase separation tends toward zero at large $N$, regardless of the variance of the random interactions. The ensemble-averaged phase behavior is presented in Figure 3a, where ${\langle\theta\rangle}$ is shown relative to the largest possible angle of phase separation for a given number of components, $\theta_{N}$. The cross-over between demixing and condensation, where $\langle\theta\rangle$ starts to tend towards 0, occurs at larger $N$ in random-mixture ensembles with greater values of $\sigma$. In Figure 3a, the cross-over occurs near ${N\simeq 4}$, ${N\simeq 16}$ and ${N\simeq 64}$ for ensembles with ${\sigma=\bar{\epsilon}/4}$, ${\sigma=\bar{\epsilon}/3}$ and ${\sigma=\bar{\epsilon}/2}$, respectively. The variance in $\theta$ across realizations of the interaction matrix also decreases as $\langle\theta\rangle$ approaches zero, while the larger variance in the ensemble-averaged angle of phase separation near ${\langle\theta\rangle\simeq\nicefrac{{1}}{{2}}}$ is a consequence of the bimodal distribution shown in Figure 2. We shall discuss the origin of this bimodality and the scaling of $N$, $\sigma$ and $\bar{\epsilon}$ at the cross-over between demixing and condensation using mean-field calculations below.

Multiphase coexistence correlates with the angle of phase separation

Free-energy landscapes calculated at the cloud point indicate that multiphase coexistence is most likely to occur in concert with demixing transitions. In Figure 3b, we plot the probability, calculated separately for each random-mixture ensemble, of finding more than one high-concentration phase in the vicinity of the cloud point. The correlation of the number of free-energy basins identified in our simulations with the ensemble-averaged angle of phase separation implies that the presence of a nearby multiphase coexistence region can be predicted by the angle of phase separation at the cloud point. This observation allows us to classify the likely phase diagram of a multicomponent mixture based on the statistical distribution of the inter-molecular interactions in a random mixture ensemble. Furthermore, we find that condensed and demixed free-energy basins only appear together in ensembles where ${\langle\theta\rangle\simeq\nicefrac{{1}}{{2}}}$. This feature indicates that the cross-over from demixing to condensation is a sharp, first-order transition that depends on the relative stability of competing high-concentration phases; the smooth transition in Figure 3a is the result of averaging over the bimodal distribution of angles of phase separation.

Mean-field theory of multicomponent

phase separation

We can gain further insight into the transition between these two opposing types of phase behaviors with a mean-field analysis of the random-mixture model. In particular, this approach reproduces the bimodality of the phase behavior observed in our simulations. This theory also allows us to predict the scaling of the condensation–demixing cross-over as a function of $N$, $\sigma$ and $\bar{\epsilon}$.

With all components present in equal concentrations in the parent phase, and ignoring spatial correlations between the various particles and vacancies on the lattice, the Helmholtz free energy of a mixture with concentration vector ${\bm{\rho}}$ and lattice coordination number $z$ is

[TABLE]

where ${\phi\equiv\sum_{i}\rho_{i}\leq 1}$ and ${\beta\equiv(k_{\text{B}}T)^{-1}}$ is the inverse temperature. In Eq. 3, the second term accounts for the mixing entropy of vacancies. Because calculating the total concentration of each phase at coexistence is challenging in the context of a multicomponent mixture, we shall make the simplifying assumption that all high-concentration phases have the same total concentration $\phi$. With this approximation, we can write the free-energy difference, ${\Delta F}$, between a phase with an arbitrary composition ${\bm{x}}$ and the equal-composition phase as

[TABLE]

where ${\Delta\epsilon_{ij}\equiv\epsilon_{ij}-\sum_{ij}\epsilon_{ij}/N^{2}}$. The first and second terms in Eq. 4 represent the differences between the two phases in the mixing entropy and the average interaction energy, respectively.

Now let us assume that the components in an arbitrary phase interact via a specific subset of $n^{2}$ interaction energies and that each of these interactions is equally likely to appear at any bond between nearest-neighbor particles. Choosing such a subset entails reshuffling the rows and columns of the interaction matrix in order to select a ${n\times n}$ sub-matrix of interaction energies, as illustrated in Figure 4a. For example, a phase with ${n=1}$ could correspond to a single-component phase, if the chosen interaction is on the diagonal of the matrix, or a two-component phase in which the components occupy alternating lattice sites, if the chosen interaction is off-diagonal. The lowest free energy of a phase with $n^{2}$ interactions, relative to the equal-composition phase with $N^{2}$ interactions, is

[TABLE]

where ${\beta^{\prime}\equiv\beta z\phi/2}$. The vectors ${\bm{\tau^{(n)}}}$ and ${\bm{\upsilon^{(n)}}\!}$, which have precisely $n$ ones and ${N-n}$ zeros, implement the reshuffling of rows and columns (Figure 4a). The second term in Eq. 5 is proportional to the maximum value (over all possible vectors ${\bm{\tau^{(n)}}}$ and ${\bm{\upsilon^{(n)}}\!}$) of the average of the $n^{2}$ interaction energies.

The composite parameter $\beta^{\prime}$ has units of inverse energy and thus sets the scale against which the standard deviation of the interaction distribution should be measured. This parameter is an implicit function of the mean interaction strength, ${\bar{\epsilon}}$, which plays an important role in determining the total concentration, $\phi$, at the cloud point. Although the functional dependence of $\phi$ on the interaction matrix is complicated and best determined through simulation (19), comparison with a one-component mean-field model indicates that $\beta^{\prime}$ is closely related to the critical interaction energy, $\epsilon_{\text{c}}$, in a single-component system. With only one component, it is easy to show that ${\epsilon_{\text{c}}=4/\beta z}$ at the critical number density ${\phi_{\text{c}}=\nicefrac{{1}}{{2}}}$, meaning that ${\beta^{\prime}(\epsilon_{\text{c}})=\epsilon_{\text{c}}^{-1}\!}$. A condensation phase transition in a multicomponent system cannot occur if ${\bar{\epsilon}<\epsilon_{\text{c}}}$. For ${\bar{\epsilon}>\epsilon_{\text{c}}}$, the single-component analogy indicates that $\beta^{\prime}$ varies slowly with ${\bar{\epsilon}-\epsilon_{\text{c}}}$. We therefore conclude that $\beta^{\prime}(\bm{\epsilon})$ in multicomponent systems should be interpreted roughly as the critical bond energy for all ${\bar{\epsilon}\geq\epsilon_{\text{c}}}$.

Mean-field predictions of condensation and demixing

Numerical analysis shows that Eq. 5 reproduces the bimodal phase behavior observed in our simulations. We evaluated Eq. 5 for many realizations of the random interaction matrix and plotted the results in Figure 4b. The second term in Eq. 5 is minimized when ${n=1}$, in which case ${\bm{\tau^{(1)}}}$ and ${\bm{\upsilon^{(1)}}}$ pick out the largest entry in the interaction matrix. As more interaction energies are incorporated, the average of the ${n^{2}}$ interactions regresses toward the mean, ${\sum_{ij}\Delta\epsilon_{ij}/N^{2}=0}$. However, this reduction in the average interaction strength is not perfectly balanced by the increase in the mixing entropy. As a result, Figure 4b shows that the free-energy difference as a function of $n$ is bistable: either a demixed phase with ${n\simeq 1}$ or a condensed phase with ${n\simeq N}$ is thermodynamically favored. Equal probabilities for these two scenarios, corresponding to coexistence between demixing and condensation phase transitions, can be obtained by tuning either the total number of components $N$ or the standard deviation of the interaction energies such that ${\Delta F=0}$ in the ${n\simeq 1}$ phase.

The free-energy curves plotted in Figure 4b resemble a classical first-order phase transition between stable demixed and condensed phases, where the total number of components and the standard deviation of the interaction strengths are the control parameters. In random mixtures with uncorrelated interaction energies, high-concentration phases with an intermediate number of components — corresponding to phases with intermediate values of $\theta$ — are unlikely to be observed. The emergence of a condensation phase transition, where ${n\simeq N}$, can be attributed to the greater mixing entropy of the equal-composition, high-total-concentration phase. In contrast, if ${n\simeq 1}$, the thermodynamic driving force for demixing results from the extreme values of the distribution of random interactions, which implies that only a small number of distinct demixed phases, corresponding to the most strongly interacting component subsets, are likely to be stable. We find that the fluctuations in the extreme values of $\bm{\Delta\epsilon}$ are large in comparison to the standard deviation of the distribution of interaction energies, leading to a significant variance in ${\Delta F(n)}$ across realizations of the interaction matrix. As a result, the mean-field theory predicts a smooth cross-over between demixing and condensation behavior when averaged over an ensemble of random mixtures.

We can predict the scaling of the conditions for coexistence between condensation and demixing transitions from the extreme-value statistics of the interaction matrix. When ${n=1}$, the expected value of the second term in Eq. 5 is ${-\beta^{\prime}\!\sigma(\sqrt{4.73\ln N}-1)}$. It follows that the coexistence curve between condensation and demixing is ${\sqrt{\ln N}=a\beta^{\prime}\!\sigma+\sqrt{(a\beta^{\prime}\!\sigma)^{2}-\beta^{\prime}\!\sigma}}$, where ${a\simeq 1.09}$. In the large $N$ limit, the resulting scaling relationship ${\sqrt{\ln N}\sim\beta^{\prime}\!\sigma}$ is considerably steeper than the central-limit-theorem scaling ${\sqrt{N}\sim\sigma/\bar{\epsilon}}$, which was previously predicted on the basis of a stability analysis of the homogeneous parent phase (17). Consequently, the $N$–${\beta^{\prime}\!\sigma}$ phase diagram presented in Figure 4c suggests that, in mixtures characterized by large values of ${\beta^{\prime}\!\sigma}$, the number of components required to promote a condensation phase transition may be extremely large. The simulation results shown in Figure 3 are also incompatible with ${\sqrt{N}\sim\sigma/\bar{\epsilon}}$ scaling, although the available data are insufficient to confirm our mean-field prediction.

This mean-field theory can be easily extended to predict the phase behavior of multicomponent mixtures under more general assumptions regarding the inter-molecular interactions and the parent-phase composition. We have verified that the predicted bistability is robust with respect to weak correlations among the inter-molecular interactions. In particular, we checked this result for mixtures where the interaction matrix follows the Lorentz–Berthelot mixing rule (24) and thus has the form ${\epsilon_{ij}\sim\sqrt{s_{i}s_{j}}}$, where the ${\{s_{i}\}}$ are independent Gaussian random variables with the constraint ${s_{i}\geq 0}$, and the largest value of ${\epsilon_{ij}}$ is always found on the diagonal of the interaction matrix. In cases where the interaction distribution is bimodal, our predictions can be applied to the strongly interacting components, which may condense or demix independently from the other weakly interacting components. In mixtures with unequal component concentrations, Eqs. 3 and 4 must be amended to account for the decrease in entropy relative to the parent phase, ${\sum_{i}x_{i}\ln x_{i}^{(0)}\!}$. The free-energy difference for each demixed phase will then depend on the parent composition as well as the selected subset of inter-molecular interactions.

Discussion

For phase separation to be a viable mechanism of spatial organization, it must be possible to tune the stability of multiple coexisting phases via a limited number of local parameters. In addition, the phase behavior of a mixture with many components should be insensitive to fluctuations in the interactions between each and every pair of components. We have shown that a model in which the inter-molecular interactions are chosen randomly produces phase diagrams that satisfy this requirement under rather general conditions. The most probable phase behavior of such a mixture falls into one of two classes, depending on the number of components and the statistical distribution of the inter-molecular interactions. Under the conditions where a homogeneous mixture first becomes unstable, either a few components will demix, or all components will condense together. The former type of phase behavior typically supports coexistence among multiple compositionally distinct phases, given small adjustments to the component chemical potentials. This means that, when demixing occurs in a random mixture, the conditions for multiphase coexistence are both readily accessible and easily tunable.

Our results have a number of important implications for biological mixtures. One consequence of bimodal phase behavior is that only a small number of interactions need to be tuned to achieve coexistence among multiple demixed phases. In the absence of strong correlations among the inter-molecular interactions, our mean-field theory predicts that each demixed phase is unlikely to be enriched in more than a few components. The stabilities of these phases can therefore be tuned by optimizing a small number of component chemical potentials. To manipulate the phase behavior within the lifetime of a single cell, such control is likely to be achieved through the regulation of macromolecular abundances, while on evolutionary timescales, the phase behavior may be tuned via local changes to the inter-molecular interactions. In an intracellular mixture, multiphase demixing results in the formation of compositionally inhomogeneous droplets, such as spontaneous nucleolar compartmentalization (11). However, the spatial ordering of the phases within a droplet depends on the interfacial free energies, which we have not examined here.

Functional multiphase coexistence furthermore requires that demixing transitions are not hidden by a dominating condensation instability. Our mean-field theory predicts that the required number of components for condensation grows rapidly with the variance of the interaction-energy distribution. This scaling implies that demixing transitions can be observed in mixtures that contain a few thousand distinct components with an interaction-distribution standard deviation on the order of a few ${k_{\text{B}}T}$. In our simulations, we have chosen the variance of the interaction strengths in order to explore the behavior of the model over a computationally tractable number of components. However, these parameters could instead be obtained from high-throughput protein–protein interaction assays. A previous analysis of Yeast 2-Hybrid experiments (25) found that the proteome-wide distribution of nonspecific interactions is indeed quite broad, with an estimated mean and standard deviation of $-4k_{\text{B}}T$ and $2.5k_{\text{B}}T$, respectively.

Our model also shows that multiphase coexistence is not simply a consequence of having a large number of components. In contrast to the complexity that might be expected from the Gibbs phase rule, phase separation into two phases with similar compositions is a common outcome for mixtures with heterogeneous interactions. This observation helps to explain the condensation-like phase behavior that has been observed for components embedded in lipid membranes (1). Although such simple phase behavior might look like that of a single-component solution, this is not a result of the interacting components being indistinguishable. On the contrary, this behavior is a consequence of the increased mixing entropy that stabilizes homogeneous-composition phases in mixtures with a large number of distinguishable components.

Owing to the simplicity of our model and the minimal constraints imposed on the form of the interaction matrix, our results can be applied to a wide variety of biological systems, including both two-dimensional membranes and three-dimensional fluids. However, additional information about the inter-molecular interactions is required before quantitative conclusions can be drawn for any specific system (26). Whereas we have assumed uncorrelated pairwise interactions in our simulations, some classes of biomolecules, such as intrinsically disordered proteins, are more likely to interact promiscuously. In particular, experiments have implicated proteins with low-complexity sequences in some examples of in vivo phase separation (27). These considerations imply a more structured interaction matrix than the uncorrelated Gaussian ensemble. In addition, protein expression levels are typically anticorrelated with the propensity of a protein to form nonspecific interactions (25; 28; 29). While our mean-field calculations indicate that bimodal phase behavior is robust to physically realistic correlations among the random inter-molecular interactions, very strong correlations may affect this result. Further investigation is warranted to examine the effects of correlated interactions and unequal concentrations on the phase behavior of multicomponent mixtures.

In summary, we have described the phase behavior of a minimal model of a biological mixture with many components. We have shown that the general features of the phase diagrams of such mixtures depend on a small number of parameters that describe the distribution of pairwise interactions between components. With a sufficiently broad distribution of weakly correlated inter-molecular interactions, these mixtures allow for tunable coexistence among multiple phases. Our results therefore suggest that multicomponent mixtures with random interactions exhibit the essential thermodynamic features that are required for phase separation to play a functional role in a biological system. While the study of high-dimensional phase diagrams has long been a topic of theoretical interest (30; 31), recent explorations of the organization of complex biological systems — which are truly multicomponent mixtures — have brought new relevance to this problem. We hope that the general principles developed here will guide future studies in this emerging field.

Acknowledgements.

W.M.J. acknowledges support from the Gates Cambridge Trust and the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1143678.