Charge Pattern Matching as a"Fuzzy"Mode of Molecular Recognition for the Functional Phase Separations of Intrinsically Disordered Proteins

Biologically functional liquid-liquid phase separation of intrinsically disordered proteins (IDPs) is driven by interactions encoded by their amino acid sequences. Little is currently known about the molecular recognition mechanisms for distributing different IDP sequences into various cellular membraneless compartments. Pertinent physics was addressed recently by applying random-phase-approximation (RPA) polymer theory to electrostatics, which is a major energetic component governing IDP phase properties. RPA accounts for charge patterns and thus has advantages over Flory-Huggins and Overbeek-Voorn mean-field theories. To make progress toward deciphering the phase behaviors of multiple IDP sequences, the RPA formulation for one IDP species plus solvent is hereby extended to treat polyampholyte solutions containing two IDP species. The new formulation generally allows for binary coexistence of two phases, each containing a different set of volume fractions $(\phi_1,\phi_2)$ for the two different IDP sequences. The asymmetry between the two predicted coexisting phases with regard to their $\phi_1/\phi_2$ ratios for the two sequences increases with increasing mismatch between their charge patterns. This finding points to a multivalent, stochastic,"fuzzy"mode of molecular recognition that helps populate various IDP sequences differentially into separate phase compartments. An intuitive illustration of this trend is provided by Flory-Huggins models, whereby a hypothetical case of ternary coexistence is also explored. Augmentations of the present RPA theory with a relative permittivity $\epsilon_{\rm r}(\phi)$ that depends on IDP volume fraction $\phi=\phi_1+\phi_2$ lead to higher propensities to phase separate, in line with the case with one IDP species we studied previously. ...

with a relative permittivity r (φ) that depends on IDP volume fraction φ = φ 1 +φ 2 lead to higher propensities to phase separate, in line with the case with one IDP species we studied previously. Notably, the cooperative, phase-separation-enhancing effects predicted by the prescriptions for r (φ) we deem physically plausible are much more prominent than that entailed by common effective medium approximations based on Maxwell Garnett and Bruggeman mixing formulas. Ramifications of our findings on further theoretical development for IDP phase separation are discussed.

Introduction
Nearly two decades of increasingly intensive research established that intrinsically disordered proteins/protein regions (abbreviated collectively as IDPs here) serve many important biological functions, and are especially critical for signaling and regulation in multicelluar organisms [1][2][3][4][5][6][7][8][9][10][11][12]. Recently, it was discovered that IDPs function not only at the level of individual molecules. An expanding repertoire of IDPs have been seen to undergo liquid-liquid phase separation in vitro, intriguingly parallelling the formation of many types of condensed liquid/gel-like bodies/organizations in living organisms, including extracellular materials, transcription complexes, nucleating sites of intermediate filament organization, and various membraneless organelles. It is apparent from these experimental observations that IDP condensation constitutes a major physical underpinning of these condensed bodies, which serve as hubs for specifically regulated sets of biomolecules to interact. As such, IDP phase separation is one of Nature's means to achieve the spatial and temporal compartmentalization necessary for the organization of vital processes [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. Although much detail remains to be ascertained, examples of membraneless organelles and an IDP species whose phase separation has been found to be a likely contributor to their assembly include chromatoid bodies, nuage or germ granules in mammalian male germ cells [33] and the DEAD-box RNA helicase Ddx4 [20], the Caenorhabditis elegans germline-specific perinuclear RNA granules known as P-granules [34] and the Ddx3 RNA helicase LAF-1 [21], as well as the stress granules triggered by integrated stress response [35] and the RNA-binding heterogeneous nuclear ribonucleoprotein A1 (hnRNPA1) [22]. Because of the importance of these bodies to biological regulation, malfunctioning of the corresponding IDP phase separation processes can lead to deregulation and diseases, including cancer due to loss of regulation of stress granules [36], protein fibrillization and thus amyloid diseases [22], and various forms of neurological disorder [37][38][39].

Seeking "sequence-phase" relationships
With the advent of IDPs, the molecular biology paradigm of seeking "sequencestructure" relationships for globular proteins has to be expanded to encompass "sequence-ensemble" relationships for individual IDPs [40][41][42]. Now, the additional question we need to ask is: How do the phase behaviors of IDPs depend on their amino acid sequences? In other words, what are the "sequence-phase" relationships [43]? Although computational study of IDPs is still in its infancy, much insight into the conformational properties and binding energetics of individual IDPs has been gained by explicit-chain simulations [8,9,[44][45][46][47][48][49]. In contrast, because IDP phase separation is a multiple-chain property, computationally it is extremely costly to simulate using fully atomic explicit-chain models [50], notwithstanding promising progress made by coarsegrained approaches that treat groups of amino acid residues of IDPs as interaction modules in continuum space [51] or on lattices [27] and simulation algorithms developed recently [52,53] for phase separations of globular proteins [54]. In this context, analytical theories of IDP phase separation are valuable not only because of their tractability, but also-and more importantly-for the conceptual framework they offer for understanding a highly complex phenomenon.

Mean-field and random-phase-approximation (RPA) theories of phase coexistence
Developed mainly for synthetic polymers at its inception, the basic statistical mechanical framework of Flory-Huggins (FH) theory [55,56] is useful for describing phase separation in the biomolecular context [17,20,57,58]. FH assumes that the interactions among the monomers (residues) of the polymers have short spatial ranges. In view of the fact that not only short-range interactions such as solvent-mediated hydrophobic effects but long-range Coulomb forces are important for driving the phase separation of certain IDPs, it has been suggested [58] that the Overbeek-Voorn (OV) theory [59,60] should be more appropriate in those cases. Inasmuch as sequence dependence is concerned, however, both FH and OV are mean-field theories that account only for composition but not sequence information. In these theories, all residues belonging to any given set of chain sequences are allowed to interact on equal footing irrespective of the correlation dictated by chain connectivity. Therefore, to address sequence specificity of IDP phase separation, one needs to go beyond FH and OV. Accordingly, we recently put forth a random-phase approximation (RPA) theory that approximately accounts for the effects of arraying different charge patterns along the IDP chain sequence [61][62][63]. The term "RPA" was first introduced by Bohm and Pines in their quantum mechanical collective description of electron interactions [64]. The analogy of this approach with the approximate polymer theory that considers terms up to quadratic in particle density was recognized by de Gennes [65]. As detailed elsewhere [62], the RPA formulation we developed [61], which follows largely that of Olvera de la Cruz [66,67], is successful in providing a physical rationalization [61,62] for the experimental salt dependence of Ddx4 phase separation as well as the difference in phase behavior between the wildtype and a charge-scrambled variant of Ddx4 [20]. Our RPA theory suggests further that the tendency for a collection of IDP chains with a given charge sequence to phase separate is strongly-but negatively-correlated with the conformational dimensions of individual IDP molecules of the sequence [63], and that both of these properties are well correlated with the charge pattern parameter κ of Das and Pappu [68] and the "sequence charge decoration" parameter SCD of Sawle and Ghosh [69]. In principle, these predictions are now testable using experimental techniques similar to those employed to study "IDP polymers" [70].

Sequence-dependent multiple-component IDP phase separation
Membraneless organelles are complex functional units consisting of many protein and nucleic acid components [18]. Different types of such units are enriched with different varieties of proteins and nucleic acids. Some individual membraneless organelles have mesoscopic substructures with different degrees of fluidity, as in the case of stress granules [71]. In a similar vein, immiscible liquid phases have been shown to contribute to subcompartmentalization of the nucleolus [27]. Clearly, a viable spatial organization of cellular processes necessitates a heterogeneous distribution of different biomolecular components into different membraneless organelles and their substructures, rather than having all IDP species condensing into a big gemisch. How is this achieved physically?
One aspect of this question was addressed recently using a simple cubic lattice model of multicomponent mixtures confined to a 6×6×6 box, with each component represented by a bead on the lattice [72]. By considering hypothetical intercomponent interaction strengths (contact energies), the authors found that with sufficient heterogeneity in contact energies, demixed domains are likely to segregate. This and other results of this big-picture study suggest that phase separation into cellular compartments with different compositions is a robust consequence of interaction heterogeneity among biomolecules [72,73].
With this in mind, the logical next step in our pursuit of sequence-phase relationships is to ascertain how interaction heterogeneity is encoded genetically. For globular protein folding, the fact that the amino acid alphabet is finite [74] implies that there are physical limits to interaction heterogeneity and structural encodability, as has been illustrated by simple exact models [75,76]. Similarly, physical limits should exist in the ability of different IDP sequences to demix. Taking a step toward deciphering what is physically achievable, here we present an RPA formulation for the phase behavior of two charged sequences as models for two IDP species. Consistent with physical intuition, we found that the tendency for the two IDP species to demix in two coexisting phases increases with increasing difference in their charge patterns. This phenomenon represents a statistical, multivalent mode of molecular recognition for cellular organization that differ from the structurally highly specific form of recognition among folded proteins but share similarities with the "fuzzy complexes" [77][78][79] involving individual IDP molecules [80][81][82][83].
Delving deeper into the role of electrostatics in IDP phase separation, we have also extended the two-sequence RPA formulation to address how a relative permittivity r (φ) that depends on IDP volume fraction φ may affect IDP phase properties. Several common effective medium approximations [84] posit a gradual decrease in r from the r (φ = 0) ≈ 80 value for pure water with increasing φ. But physical consideration [62] and experimental volumetric measurements suggest a much sharper decrease, with r (φ = 0.2) ≈ 20, leading to large cooperative effects that enhance phase separation significantly. These findings and their ramifications are detailed below.

Theoretical development of the RPA formulation
The development of RPA theory for a pair of charged sequences constitutes the bulk of the results presented in subsequent sections of this article. This effort is based on an extension of the RPA formulation for a single sequence that we put forth recently [61,62].

Experimental determination of dissolved protein volumes
To address the effect of volume of dissolved proteins on the relative permittivity of the resulting aqueous solution, nuclear magnetic resonance (NMR) and absorbance measurements were performed on two folded globular proteins bovine serum albumin (BSA) and hen egg white lysozyme (HEWL) as well as two IDPs Ddx4 14FtoA and Ddx4 cond , which are, respectively, a mutant of Ddx4 in which all 14 phenylalanines are mutated to alanines and the concentrated phase of phase-separated wildtype Ddx4.

Measurement of water content of protein samples
NMR spectra were recorded using a Bruker Ascend III spectrometer at 14.0 T equipped with a cryogenically cooled triple resonance gradient probe. Spectra were processed using NMRPipe [85]. 1D 1 H spectra were recorded on protein samples over a range of concentrations (5-400 mg mL −1 ) and the integrated water signals were compared with the corresponding integrals obtained from a spectrum recorded of buffer (the same buffer composition as used for each protein sample). BSA was purchased from Sigma and samples were prepared in 20 mM sodium phosphate (NaPi), 100 mM sodium chloride (NaCl), 10 % 2 H 2 O/90 % 1 H 2 O, pH 6.5. HEWL was purchased from BioBasic and samples were dissolved in 20 mM sodium citrate, 100 mM NaCl, 10 % 2 H 2 O/90 % 1 H 2 O, pH 5 (lower pH was used due to limited solubility of HEWL in NaPi at pH 6.5). Ddx4 14FtoA samples were prepared according to [86] and dialysed against 20 mM NaPi, 100 mM NaCl, 5 mM tris(2-carboxyethyl)phosphine (TCEP), 10 % 2 H 2 O/90 % 1 H 2 O, pH 6.5. For Ddx4 cond the same buffer was used but the NaCl concentration was varied between 100-400 mM in order to generate samples with protein concentrations between 200 and 400 mg mL −1 . For phase-separated samples, it was ensured that the entirety of the probe coil was occupied by the condensed phase, thus avoiding contaminating signals from the more hydrated dilute phase. Spectra were recorded using both small flip angle (θ < 10 • ) and θ = 90 • pulses with very similar results in both cases.

Overview of three-component phase behaviors
Possible phase behaviors of a three-component liquid system are outlined in Fig. 1. In general, the system can be a homogeneous solution [ Fig. 1(b], or it can separate into two coexisting phases [binary coexistence; Fig. 1(c)-(e)], or separate into three coexisting phases [ternary coexistence; Fig. 1(f)]. Fundamentally, phase behavior is governed by the intra-and inter-component interactions as well as environmental conditions such as temperature and pressure. Our theories below provide a rudimentary physical account of how interactions among IDP chains with two different amino acid sequences affect the conditions under which binary and ternary coexistence emerge. The theories presented here are for solution systems with an effective infinite volume. As such, our theories account for the differences among scenarios typified by the leftmost drawings in Fig. 1(c)-(f) but they are not equipped to address details such as droplet size and geometry. In other words, they provide no discrimination among different droplet geometries along a given horizontal row in Fig. 1. Accounting for the latter would require additional modeling of the interfacial tensions between different solution phases [27].

RPA theory for two charged sequences
Based on our previous RPA formulation for a single sequence [61,62], the approach is now extended to consider two model polypeptide (IDP) sequences s 1 and s 2 in a salt-free aqueous solution. Using notation similar to before [62], the electric charges along the sequences are written as 2 , . . . , σ (N 2 ) 2 }, where N 1 and N 2 are the numbers of residues in s 1 and s 2 , respectively. The corresponding volume fractions of the IDPs in solution are denoted as φ 1 and φ 2 . The present study is restricted to IDP sequences with zero net charge. Counterions are not considered.

Free energy as a function of two IDP volume fractions
Following the FH lattice argument [55], we partition the spatial volume V of the solution system into lattice units, a 3 , that corresponds to the volume of a solvent molecule. Accordingly, the RPA free energy per unit volume and in units of the product of Boltzmann constant k B and absolute temperature T is cast as the per-lattice-site quantity . (e) One phase is concentrated in seq2 but dilute in seq1 (top), the other phase is concentrated in seq1 but dilute in seq2 (bottom). (f) Ternary coexistence. In this example, the system separates into a dilute phase for both seq1 and seq2 (top), a phase concentrated in seq2 but dilute in seq1 (middle), and a phase concentrated in seq1 but dilute in seq2 (bottom), with water present in all three phases. Circle(s)-in-square drawings to the right of (c)-(f) depict possible configurations of phase-separated droplets in the given scenario. In each case, the rightmost schematic phase diagram illustrates the manner in which a given set of bulk concentrations (bulk volume fractions for seq1 and seq2) of (φ 1 , φ 2 ) represented by the central dot with the color of (b) is separated (arrows) into two or three phases with different (φ 1 , φ 2 )'s. The volume fraction of water is equal to where the negative entropy −s is the entropic contribution to free energy in units of k B T . This term is given by the standard FH entropy of mixing for a system comprising of s 1 , s 2 , and solvent: where 1 − φ 1 − φ 2 is the volume fraction of solvent. The electrostatic contribution f el is calculated by RPA [61,62], viz., where k is the reduced wave number that absorbs the virtual bond length b a of the polypeptide backbone by re-defining kb =k in Ref. [62] as k (i.e., kb → k), such that where 4π/[k 2 (1 + k 2 )] is from the Fourier transformation of Coulomb interaction with a short-range cutoff [61,67] in units of k B T , 0 is vacuum permittivity and r is relative permittivity, T * ≡ b/l B is the reduced temperature defined by Bjerrum length l B = e 2 /(4π 0 r k B T ). Here |q is the (N 1 + N 2 )dimensional column vector representing the two charge sequences, namely q i = σ for N 1 + 1 ≤ i ≤ N 1 + N 2 , q| is the transposed row vector, q|Ĝ k |q ≡ ij q i (Ĝ k ) ij q j with (Ĝ k ) ij being the i, j element of the bare two-body correlation matrixĜ k of all possible sequence-sequence correlations [62], andĜ 12 (k) =Ĝ 21 (k). As in our previous studies [61,62], we consider a simple formulation in which all IDPs are modeled as Gaussian chains without excluded volume within the RPA formalism [66]. In this approximation, there is no correlation between the positions of different chains; henceĜ 12 (k) =Ĝ 21 (k) = 0 and follow from the average of exp(ik · R ij ) over a Gaussian chain ensemble wherein R ij is the vector between chain positions i, j and k 2 = k · k (Eq. (IX.59) of [65]). Eq. (4) then becomes where

Free energy landscape and spinodal instability
As an example, we first apply this formulation to two N 1 = N 2 = 50 charged IDP sequences corresponding to sv28 and sv24 in Das and Pappu [68]. The sequences are labeled here as seq1 and seq2 respectively. The RPA free energy function for the two sequences ( Fig. 2) consists of regions of different curvatures: parts of the landscape are convex downward (i.e., has a convex downward curvature) whereas some other parts are convex upward (concave downward). If a given IDP solution has overall (bulk-state) IDP volume fractions (φ 0 1 , φ 0 2 ) situated in a convex-downward region, the state is thermodynamically stable and thus the bulk-state volume fractions are maintained. In contrast, if the bulk state is in a concave-downward region, the state is thermodynamically unstable because a phaseseparated state allows for a lower free energy. Accordingly, the system undergoes phase separation to multiple coexisting phases with different protein volume fractions [62].
The boundary between the convex and concave regions is determined by the saddle point condition [62,88] detF ≡ This boundary defines a spinodal instability region in which detF < 0. Because the determinant of a matrix is equal to the product of all its eigenvalues, this instability condition indicates that one, but not both, of the eigenvalues ofF is negative. This means that second-order perturbations of free energy with respect to φ 1 , φ 2 along the direction of the corresponding eigenvector diverges, signaling that the system cannot maintain a homogeneous phase. Note that although the opposite of the detF < 0 condition, viz., detF > 0, does not by itself exclude the possibility that both eigenvalues ofF are negative and thus the system is thermodynamically unstable despite not satisfying detF < 0 (e.g. at a local maximum), exhaustive numerical searches (φ 1 , φ 2 = 0.001, 0.002, . . . , 0.999) did not find any such instance for all RPA and FH systems studied in this work. Hence, for these systems, Eq. (10) is the valid condition for spinodal boundaries, examples of which are shown as dashed curves in Fig. 2 for the (seq1, seq2) system. shown along the axes in (a) and (c), wherein positive and negative charges are depicted as red and blue beads respectively, correspond to sv28 and sv24 in Das and Pappu [68]. φ 1 and φ 2 are their volume fractions, respectively, as the sequences are re-labeled as seq1 and seq2 in this work (Table 1). Results in this figure are for T * = 4. (a) Free energy landscape. The plotted quantity is the free energy f (φ 1 , φ 2 ) in Eq. (1) minus a linear function a 1 φ 1 + a 2 φ 2 of φ 1 , φ 2 where the coefficients are chosen to be a 1 = −1.1447 and a 2 = −1.1453. The sole purpose of subtracting this linear term is to graphically highlight the changes in curvature on the landscape. The subtraction has no effect on the determination of phase coexistence (Sec. 4.3). As examples of phase coexistence, three pairs of phase-separated states are marked by grey dots, whereby each pair is connected by its own tie line (grey solid line). (b) View of the landscape in (a) from an elevated vantage. The two dashed curves are spinodal phase boundaries defined by Eq. (10) (one of the boundaries is very close to the origin with φ 1 , φ 2 intercepts ≈ 0.003, see zoom-in view in (d)); the region between these two boundaries satisfies detF < 0. The three pairs of dots marking coexisting phases connected by tie lines are the same as those in (a). Note that subtracting a linear function of φ's from f (φ 1 , φ 2 ) does not alter the spinodal boundaries determined by the matrixF that contains only second derivatives of φ's. (c) Heat map representation (scale on the right) of the function 10 4 × [f (φ 1 , φ 2 ) − a 1 φ 1 − a 2 φ 2 ] in (a). Pairs of coexisting (φ 1 , φ 2 ) connected by tie lines and spinodal boundaries are the same as those in (b). (d) Zoom-in plot of the 0 ≤ φ 1 , φ 2 ≤ 0.005 region in (c) to provide a clearer view of the spinodal phase boundary that is very near to the origin (dashed line) and the dilute phases of the three coexisting pairs.
In two-component systems such as those consisting of a single IDP species and solvent molecules, spinodal instability necessarily leads to binary coexistence [62]. In systems with more than two components, more coexisting phases are allowed. According to the Gibbs phase rule, the maximum number of coexisting phases under a given set of environmental conditions is equal to the number of components in the system [62,89] but, depending on system specifics, the actual number of coexisting phases can be smaller. Spinodal instability always implies that the system existing as a single phase is thermodynamically untenable and thus phase separation must occur. However, when the number of components is larger than two, the precise number of coexisting phases is governed by how total free energy varies with changes in the component volume fractions in different combination of phases or, equivalently, by the balance of chemical potentials for each component across different phases [62,88]. As examples, three instances of binary coexistence in the (seq1, seq2) system are depicted in Fig. 2 as pairs of dots connected by solid tie lines. The general procedure for determining the conditions for such coexistence is as follows.

Binary coexistence in three-component systems: Applications to IDPs with two different sequences plus solvent
As mentioned, we use two methods to determine phase equilibrium: by ascertaining the minimum free energy among single-and multiple-phase states [90], and by balancing the chemical potentials for each of the components across different phases [56,59,91]. The two approaches are mathematically equivalent. They can be applied simultaneously to yield more accurate numerical results. Here we describe in detail how the two approaches are applied to study three-component systems that undergo binary phase separation.
A system is dictated by thermodynamics to seek its lowest free energy state. Whether a system phase separates can be ascertained by comparing its overall free energy with and without phase separation. For a system with bulk IDP volume fractions (φ 0 1 , φ 0 2 ), the free energy f bulk without phase separation and the free energy f sep for separating into two phases (labeled as α and β) that take up fractional volumes v α = v (0 ≤ v ≤ 1) and v β = 1 − v of the total system volume and with IDP volume fractions , respectively, are given by wherein conservation of volume of each of the components implies that By rewriting Eq. (12) to express IDP volume fractions in β as functions of v with 0 < v < 1 and volume fractions in α: Note To find the set of variables that yields the global minimum of f sep , we numerically search the three-dimensional space of (v, φ α 1 , φ α 2 ) by implementing the sequential least squares programming (SLSQP) algorithm [92] using the scipy.optimize.minimize function in Scipy, a Python-based numerical package for scientific computation [93]. If a given , the system is judged to be in a state of binary phase separation to the two phases α and β. In contrast, if all f sep for a given (φ 0 1 , φ 0 2 ) are larger than f bulk , the bulk state (φ 0 1 , φ 0 2 ) is thermodynamically stable and the system does not phase separate. Unlike in the two-component case in which the two separated phases are unique and independent of the bulk IDP concentration/volume fraction insofar as it is in the phase-separated regime, in three-component systems different bulk IDP concentrations/volume fractions can result in different α, β phases. Thus a complete binary phase diagram is generated by considering all possible (φ 0 1 , φ 0 2 ) combinations. For two-component systems (one IDP sequence plus solvent), we have shown that linear terms of φ in the system free energy do not affect the determination of phase equilibrium [62]. In the same vein, here we demonstrate that the same principle applies also to three-component (two IDP sequences plus solvent) systems. As described above, binary coexistence is governed by the free energy difference Consider a modified free energy g(φ 1 , φ 2 ) with an additional arbitrary linear function a 0 + a 1 φ 1 + a 2 φ 2 of φ's where a 0 , a 1 , and a 2 are constants: Now, the free energy difference between phase-separated and bulk phases becomes because a 0 in g sep and g bulk cancel. The two bracketed terms in Eq. (17) Table 1. Sequences studied are identified as seq1-7 in this article. They correspond to seven of the thirty 50-residue charged sequences with zero net charge in Das and Pappu [68]. The sequences' sv labels, the values of their charge pattern parameter κ and simulated single-chain radius of gyration R g are those in the same reference [68]. The charge pattern parameter SCD is that of Sawle and Ghosh [69]. Values of SCD and the RPA-predicted critical temperature T * cr are from Lin and Chan [63].
added to f , such as the one utilized in Fig. 2 for graphical clarity, has no impact on phase separation. We apply the above-described minimization procedure for three-component (two IDP sequences plus solvent) systems to determine (α, β) for selected pairs of sequences in Table 1. To minimize possible numerical errors, every set of {φ α 1 , φ α 2 , φ β 1 , φ β 2 } obtained by minimizing f sep is subject to further testing by comparing the chemical potentials in α and β. As described in Eq. (A.5) of Ref. [62], phase equilibrium implies the following equalities, and is the chemical potential of water [62]. Making use of the volume conservation conditions in Eq. (18) and substituting Eq. (13a) for φ β 1 and Eq. (13b) for φ β 2 , Eq. (18) becomes three equalities for three variables φ α 1 , φ α 2 , and v. It follows that a unique determination of the phase-separated volume fractions φ α 1 , φ α 2 , φ β 1 , and φ β 2 is afforded by Eq. (18). It is straightforward to show that the set of phase-separated volume fractions (18) are identical to that obtained by minimizing f sep in Eq. (14). A necessary condition for the minimization of f sep is that its Jacobian vector J sep of first-order partial derivatives of independent variables vanishes: In other words, wherein we have utilized Eq. (13a) for φ β 1 and Eq. (13b) for φ β 2 . Clearly, Eqs. (22a) and (22b) are equivalent to Eqs. (18a) and (18b), respectively, and Eq (22c) is equivalent to Eq. (18c) by virtue of Eq. (20). Q.E.D.
Starting with (α, β) obtained by minimizing f sep in Eq. (14), only those that deviate less than 0.1% from the chemical-potential-balancing equalities in Eq. (18) are accepted as valid binary pairs in our analysis. For the sequence pairs (seq1, seq5) and (seq1, seq6), a smaller threshold of 0.01% is used to ensure accuracy of the computed phase-separated φ's because for these sequence pairs the chemical potential balancing conditions are quite insensitive to variations of the φ's.

Binary coexistence of two charged sequences
Using RPA, we investigated previously how the phase separation behaviors of charged IDP sequences are affected by their charge patterns [61][62][63]. In particular, for the set of thirty KE sequences of Das and Pappu [68] with zero net charge but an equal number of 25 positively charged lysine (K) and 25 negatively charged aspartic acids (E) in different permutations, the critical temperature T * cr of phase separation was found [63] to be correlated with charge pattern parameters κ [68] and where i, j label the residues with charges σ i , σ j along a chain of length N [69]. The κ and SCD parameters exhibit similar correlations with single-chain radius of gyration R g [68,69]. The correlation of T * cr and R g with SCD is stronger than that with κ. A likely reason is that SCD accounts for nonlocal effects between charges far apart along the chain sequence whereas κ does not [63]. Although we use only SCD in our analysis below, an equivalent analysis using κ is expected to produce a similar trend.
How does the phase behavior of an IDP solution with two sequences depend on the sequences' difference in charge patterns? Intuitively, when two sequences with different SCD values are present together, their different propensities to phase separate are expected to interfere. Indeed, such an effect of inter-sequence interference is seen clearly in the Taylor expansion of the integrand of the RPA expression f el for the electrostatic contribution to free energy in Eq. (3), whereĜ * 11 (k) andĜ * 22 (k) are the product of 4π/[k 2 (1 + k 2 )T * ] with, respectively, thê G 11 (k) andĜ 22 (k) in Eq. (8). The first two terms after the second equality in Eq. (24) are self-interactions of the two sequences, identical to those in one-sequence RPA theory (see, e.g. Eq. (1) in Ref. [63]). The third term represents the interference effect in RPA. Since it is the product of square roots of the two self-interaction terms, its strength is intermediate between them, suggesting that phase behaviors of two-sequence systems are sensitive to the similarity/dissimilarity in charge pattern between the two sequences.
To investigate this sensitivity, we use the sequences in Table 1 to compute the phase diagrams of six pairs of sequences, namely seq1 with each of the six other sequences. The pairs are selected to represent a broad range of similarity/dissimilarity in charge pattern as quantified by the difference in SCD values: from the (seq1, seq2) pair with SCD = (−15.99, −17.00) to (seq1, seq7) with SCD = (−15.99, −0.41). To compare the phase behaviors of the six sequence pairs on an equal footing, all phase diagrams in Fig. 3 are computed at the same reduced temperature T * = 4. Noting that condensedphase volume fractions tend to decrease with increasing T * , this temperature is chosen because it falls in the mid-range of the broad span of T * cr 's for the sequences in Table 1. T * = 4 is much higher than the T * = 0.55 equivalent of room temperature (T = 300 K) when an aqueous r = 80 is assumed [63]. This seemingly unphysical condition in our calculation has little impact, however, on the present goal of ascertaining general principles and behavioral trends. Although we do not aim for direct, detailed comparison with experiment here, T * = 4 is experimentally relevant, for example, to IDPs with charge patterns similar to those considered here but with their electrostatic interaction strength significantly scaled down for various physical reasons such as screening or a more sparse charge distribution along the IDP sequence.
For each of the six phase diagrams in Fig. 3, the area surrounded by blue dots and "shaded" by black dashed lines is the region of binary coexistence. In other words, bulkstate volume fractions falling within this region will phase separate into two coexisting liquid phases [as in Fig. 1(c)-(e)], whereas bulk-state volume fractions residing outside this region will be stable as a single liquid phase [as in Fig. 1(b)]. Every black dashed  (c)    Table 1. Horizontal axes (φ 1 ) refer to the volume fraction of seq1, vertical axes (φ 2 ) are for the volume fraction of the other sequence in each of the pairs. The charge patterns of the sequence pairs vary from being very similar (a) to very dissimilar (f). This trend is quantified by the difference in −SCD values between the two sequences in each of the six pairs. Each dashed line is a tie line connecting a pair of blue dots that represent coexisting phases α = (φ α 1 , φ α 2 ) and β = (φ β 1 , φ β 2 ). All bulk volume fractions (φ 0 1 , φ 0 2 ) lying on one tie line undergo phase separation to the same (α, β); and, following Eq.
, which is equal to the length ratio of the two segments of the tie line from the bulk volume fractions (φ 0 1 , φ 0 2 ) on the phase diagram to the phase boundary marked by the blue dots. The black inclined solid line in (f), φ 1 + φ 2 = 1, delimits the region φ 1 + φ 2 ≤ 1 within which IDP volume fractions may vary. Insets in (a) and (b) are zoom-in plots of a part of the phase diagram with extremely low φ 1 and φ 2 . They offer a clearer view of the dilute phase boundaries for the (seq1, seq2) and (seq1, seq3) pairs. Insets in (e) and (f) are zoom-in plots of the grey-shaded regions of the respective phase diagrams. To provide a scale for comparison, the size of the grey-shaded regions φ 1 , φ 2 ∈ [0, 0.06] in (e) and (f) is chosen to be equal to the plotted regions for the other four sequence pairs. line is the tie line connecting a pair of blue dots representing separated phases (α, β) for any bulk-state volume fractions lying on the given tie line [corresponding to the arrows in the schematic phase digrams for Fig. 1(c)-(e)].
The average slope of the tie lines changes from positive for similarly patterned sequences [ Fig. 3(a), (b)] to negative for very differently patterned sequences [ Fig. 3(e), (f)]. This trend may be understood as follows. When both sequences of a given pair can undergo phase separation individually (T * cr > 4 for both), the tie lines near the φ 1 and φ 2 axes must be close to being parallel to the axes because when either φ 0 1 → 0 or φ 0 2 → 0, the two-sequence system reduces to the corresponding single-sequence system that phase separates. This situation applies to (seq1, seq2) and (seq1, seq3), resulting in positive tie-line slopes, indicating that the populations of the two sequences in each pair are well mixed even when they undergo phase separation. They prefer to stay together after they phase separate, with similar population ratios for the two sequences in the "both-dilute" (small φ's) as well as the "both-condensed" (larger φ's) phases [as in Fig. 1(c)].
Because T * cr < 4 for the other four sequences (seq4-seq7), they do not phase separate by themselves individually and therefore tie lines near the φ 2 -axis need not be approximately parallel to it. Nonetheless, tie lines close to the φ 1 -axis are still required to essentially line up with the axis. For tie lines that possess large φ 2 values, the volume conservation condition φ 1 + φ 2 = 1 enforces negative tie-line slopes. The combined effect of these constraints lead to tie-line slopes that gradually change from ≈ 0 near the φ 1 -axis to ≈ −1 near the φ 1 + φ 2 = 1 boundary, as exemplified by the case of (seq1, seq7) in Fig. 3(f). As shown in Fig. 3(d) and (e) for (seq1, seq5) and (seq1, seq6), this trend is apparent even when the phase-separated regime does not extend all the way to the φ 1 + φ 2 = 1 boundary. Negative tie-line slopes imply various degrees of demixing of the populations of the two sequences: the phase-separated state (α, β) now comprises one φ 1 -enriched (φ α 1 φ α 2 ) phase coexisting with one φ 2 -enriched (φ β 2 φ β 1 ) phase. The degree of population demixing depends on how dissimilar are the charge patterns of the two sequences in the pair. For large differences in SCD as in Fig. 3(f), one of the coexisting phases can have a very low population of seq1 but a substantial seq7 volume fraction, whereas the other phase has a relatively low population of seq7 but a substantial seq1 volume fraction [as in Fig. 1(e)].
The (seq1, seq4) pair in Fig. 3(c) is at the crossover between the well-mixed and demixed extremes. Tie-line slopes in this case are all ≈ 0, indicating that although increasing seq4 volume fraction decreases the phase separation tendency of seq1, even in the phase-separated regime the concentration of seq4 is essentially identical in the two coexisting phases, i.e., φ α 2 φ β 2 φ 0 2 [as in Fig. 1(d)]. The difference in the ratio of component populations in coexisting phases (α, β) may be quantified by comparing the volume ratio of the two sequences in the two phases.  Figure 4. Trend of the component ratio φ 1 /φ 2 in coexisting phases (α, β) for different pairs of charged sequences. Two measures, ∆ αβ (φ 1 /φ 2 ) and A αβ , are plotted. When the charge patterns of the sequences in a pair are similar (small SCD 1 −SCD 2 ), the two sequences tend to be well-mixed with similar relative volume fractions in a "both-dilute" phase and a "both-condensed" phase. When the charge patterns of the sequences in a pair are dissimilar (large SCD 1 −SCD 2 ), the two sequences tend to demix in that they largely exclude each other in the two separated phases. This mixing/demixing behavior is less extreme for intermediate SCD 1 −SCD 2 values.
We first consider a rather intuitive measure of compositional asymmetry between coexisting phases, where the absolute value ensures that ∆ αβ (φ 1 /φ 2 ) is α ↔ β symmetric, and the bracket . . . denotes averaging over all (α, β) pairs of coexisting phases. One disadvantage of this measure, however, is that the average is strongly dominated by those pairs of (α, β) with large φ 1 but small φ 2 , i.e., coexisting pairs that are close to φ 1 -axis in Fig. 3. Therefore, we also consider another composition asymmetry measure that avoids this potentially problematic feature by replacing the ratio φ 1 /φ 2 with its arctangent value normalized by π/2 such that 0 ≤ A αβ ≤ 1. Summarizing our findings using two-sequence RPA theory, Fig. 4 shows the variation of ∆ αβ (φ 1 /φ 2 ) as well as A αβ with the difference in SCD values of the six sequence pairs in Fig. 3. A reasonable correlation is seen for both composition asymmetry measures, with the A αβ measure exhibiting a better correlation by varying monotonically with SCD difference, indicating that compositional asymmetry or degree of demixing of phase-separated populations as quantified by A αβ is positively correlated with the difference in charge patterns as quantified by difference in SCD values. This plot illustrates graphically how a stochastic, multivalent form of molecular recognition that arises from the interactions among the diverse conformations in a multiple-chain ensemble can lead to demixing of different IDP species into different coexisting phases.

Comparison with phase separation in Flory-Huggins (FH) models
We next seek a deeper understanding of the RPA results and their ramifications by comparing them with the predictions of a variety of FH models. As emphasized, unlike RPA, FH by itself does not address the physics of sequence dependence [20,62]. Accordingly, FH χ interaction parameters for IDP sequences have to be provided phenomenologically by experiment or theoretically by microscopic physical theory. ‡ For example, as will be discussed further below, an intuitive and semi-quantitative connection between RPA and FH is provided by the expansion in Eq. (24). It should also be noted that FH neglects interaction terms that are higher than quadratic order in IDP volume fractions/concentrations (φ's) such as the O(G(k) 3 ) terms in Eq. (24) because G(k) ∝ φ [Eqs. (9a) and (9b)]. This approximation can be problematic when IDP concentrations are high. Nonetheless, by treating the three parameters χ 11 , χ 22 , and χ 12 in the three-component FH interaction term for two IDP species plus solvent as free (arbitrary) variables, we can either match FH behavior to that of RPA to gain conceptual insights or explore other interaction scenarios that might be physically plausible when interactions other than the rudimentary electrostatics embodied in RPA are included in the physical picture.

FH models that imitate RPA theory by having two independent χ's
We begin this analysis by first constructing FH models with interaction schemes similar to RPA, then tuning the interaction parameters to produce phase behaviors similar to those predicted by RPA in Fig. 3 T * ]} 2 such that χ 11 = dkχ 11 (k), χ 22 = dkχ 22 (k), and χ 12 = dk χ 11 (k) χ 22 (k), from which it is clear that only two set of interaction parametersχ 11 (k) andχ 22 (k) as functions of k are independent. Here we approximate this dependence by constraining χ 12 = √ χ 11 √ χ 22 . ‡ In the caption describing the FH results in Fig. 8 of Ref. [62], r d is in fact the symbol r for the equilibrium spacing in Eq.(S19) of Ref. [20]. This typographical error does not affect the results. It should also be noted that because Ref. [20] equates the ionic strength I with [NaCl] but not 2[NaCl], their effective Debye length is 4.3Å instead of the correct value of 3.04Å. To facilitate comparison with Ref. [20], however, the effective Debye length in Fig. 8 of Ref. [62] was also set to 4.3Å. We then construct three FH systems that have χ 11 = χ 22 , χ 11 χ 22 , and χ 11 χ 22 , corresponding respectively to sequence pairs with small, intermediate, and large charge pattern (SCD) differences. The phase diagrams of these models (Fig. 5) exhibit a trend similar to that seen in the RPA-predicted Fig. 3. Specifically, Fig. 5(a) is similar to Fig. 3(a), Fig. 5(b) to Fig. 3(c), and Fig. 5(c) to Fig. 3(f). This correspondence offers conceptual clarity because the degree to which the interaction between the two IDP species is favorable is explicit in FH. When χ 11 = χ 22 = χ 12 , the two species are miscible and their phase separation propensities are identical, resulting in the coexistence of one both-dilute phase and one both-condensed phase. The similarity between Fig. 5(a) Fig. 3(a) indicates that this behavior can be achieved physically by two IDP species with similar charge patterns. In contrast, when χ 11 χ 12 χ 22 , miscibility of the two species is poor and their phase separation propensities are quite different, resulting in a high degree of population demixing. The similarity of this behavior shown in Fig. 5(c) with that in Fig. 3(f) underscores once again that charge-pattern mismatches between IDPs can lead to substantially weakening of attractive interactions.

FH models with three independent χ's
We next consider the general case in which the three χ's in Eq. (27) are independent. Although this modeling setup does not have a simple correspondence with IDP sequences interacting via physical forces like that described above, the expanded variety of scenarios explored here would be valuable when behaviors much more complex than those allowed by our current RPA formulation are considered in more comprehensive and detailed physical theories. As simple examples of the rich possibilities, here we focus on FH models with χ 11 = χ 22 (≡ χ) and variable χ 12 values that are not related to χ. Fig. 6 shows the phase diagrams of three representative models with (a) χ 12 χ, (b) χ 12 χ, and (c) χ 12 < χ. Compared to the χ ≡ χ 11 = χ 22 = χ 12 case in Fig. 3(a), it is clear that a stronger inter-component attraction (larger χ 12 ) makes the phase-separated region bulge [ Fig. 6(a)], whereas a weaker inter-component attraction (smaller χ 12 ) shrinks it [ Fig. 6(b)]. Nonetheless, the tie-line slopes are positive in both situations, indicating that the two components are largely miscible. Fig. 6(b) indicates that a weakened χ 12 shrinks the phase-separated region most around φ 1 = φ 2 . As χ 12 decreases further, inter-component attraction all but vanishes, micibility disappears, resulting in the phase-separated region being broken into two parts, one for φ 1 φ 2 and the other for

Ternary coexistence in FH model
If χ 12 is made even weaker than that in Fig. 6(c), the region of poor miscibility below the φ 1 + φ 2 = 1 boundary would grow. Finally, the three phase-separated regions intersect and a new ternary coexistence region emerges in-between.
According to the Gibbs phase rule, an n-component system can separate into at most n coexisting phases, when all other environmental conditions, e.g. temperature and pressure, are kept constant [62,89]. In our system of n = 3 (two sequences plus solvent), whether the system will separate to two or three coexisting phases may be deduced by observing the variation of tie-line slopes in putative regions of binary coexistence: If the tie lines have to become parallel to the φ 1 -axis, φ 2 -axis, or the φ 1 + φ 2 = 1 line when they approach these boundaries respectively, the tie-line slopes have to be able to vary smoothly to satisfy these constraints in order for binary coexistence to be stable. In that case, ternary coexistence is unlikely. Conversely, if there are conflicts that prevent a smooth change of tie-line slope, a ternary coexistence region ensues.
An example is provided by using Fig. 6(c) as starting point. Here the three phaseseparated regions are close to the three boundaries, and their tie lines are essentially parallel to the respective boundaries. Under this circumstance, when effective intercomponent repulsion is enhanced by weakening χ 12 to cause the three regions to evolve toward merging, the conflict among the three different trends of tie-line slopes necessitates reconcilation by a region of ternary phase separation (Fig. 7).
In order to determine the three phases in ternary coexistence mathematically, we extend the phase-separated free energy expression in Eq. (11b) for phases (α, β) to including one additional phase γ, viz., where the fractional volumes v α , v β are functions of the six φ's by virtue of volume conservation, where x = 1, 2 and y is over y = α, β. Now the equalities in Eq. (18) have to include the addition phase to become Because here we have six equalities in Eq. (30) for six phase-separated φ y x 's, the solution for ternary coexistence of (α, β, γ) is unique irrespective of the bulk-state volume fractions insofar as they fall within the ternary coexistence region. This situation is different from that of binary coexistence in which the separated phases (α, β) can be different for different bulk-state (φ 0 1 , φ 0 2 )'s when they are on different tie lines. Similar to the binary coexistence case in Sec. 4.3, we proceed to demonstrate that minimizing Eq. (28) is equivalent to solving the equations in Eq. (30). Using essentially the same approach, we rewrite f ternary as a function of six independent variables: v α , v β , φ α 1 , φ α 2 , φ β 1 , and φ β 2 by first utilizing Eq. (29) to express φ γ 1 and φ γ 2 as where x, y and y have the same meanings as above. In this notation, Eq. (28) becomes Similarly to Eq. (21) for binary phase separation, we calculate the six derivatives of f ternary and set them to zero as in Eq. (22) as necessary conditions for the minimization of f ternary , resulting in where y sums over y = α, β. Substituting x = 1, 2 in Eq. (33a) yields Eqs. (30a) and (30b), whereas substituting y = α, β in Eq. (33b) yields Eq. (30c). Q.E.D. Fig. 7 provides an FH phase diagram with both binary and ternary coexistence. In this example, χ 12 is significantly smaller than χ 11 = χ 22 , resulting in strong effective repulsion between the two sequences. Consequently, the two islands of binary coexistence around the φ 1 -and φ 2 -axes intersect with the top-right binary region with (φ α 1 , φ α 2 ) = (φ β 2 , φ β 1 ), resulting in a ternary phase separation region corresponding to the αβγ triangle and its interior [marked by blue lines in Fig. 7(a) and turquoise lines in Fig. 7(b) and (c)]. The thermodynamic stability of the ternary phase-separated state within this region is illustrated by the ∆f quantity plotted in Fig. 7(b) and (c); ∆f (φ 1 , φ 2 ) is the bulk (not-phase-separated) free energy f bulk minus the free energy value for the same φ 1 , φ 2 on a plane defined by the three ternary phases (i.e., ∆f = 0 for the three points corresponding to the α, β, γ phases). Because ∆f > 0 for any other point within the triangular region, the ternary state is more stable than the bulk state in this region. Moreover, the free energy f sep of any putative binary coexistence state of a bulk state within the triangular region must lie on a tie line joining two points on the landscape in Fig. 7(b) and (c). Because ∆f > 0 for any point other than the three ternary phases in the entire plotted region-including points outside the triangular region, ∆f > 0 holds also for any putative binary coexistence state for the bulk state within the triangular region, implying that they are less stable than the ternary phase-separated state in the region.
For any given bulk-state (φ 0 1 , φ 0 2 ) in the ternary region, fractional volumes v α , v β , and v γ in the respective coexisting phases α, β, and γ are determined by solving Eq. 29 and setting v γ = 1 − v α − v β . In terms of the three-dimensional vectors , where the bracketed function linear in the φ's specifies the plane defined by α, β, and γ. In other words, the a coefficients are determined by solving f bulk (φ z 1 , φ z 2 ) = a 0 + a 1 φ z 1 + a 2 φ z 2 for z = α, β, γ. Note that ∆f > 0 for all (φ 1 , φ 2 ) in the plotted region, indicating that phase separation is preferred if the volume conservation condition for either binary [Eq. (12)] or ternary [Eq. (29)] coexistence can be satisfied. The logarithmic color scale on the right is for both (b) and (c). The turquoise lines in (b) and (c) mark the same ternary phase boundaries as those in (a). (c) Contour plot of ∆f . The three ternary coexisting phases (∆f = 0) are seen to be situated in three different basins with small ∆f .
Because the area of a triangle defined by two vectors is equal to half of the magnitude of their cross product, these fractional volumes correspond to specific ratios of triangular areas as described in the caption for Fig. 7. 6. Discussion

Insights into cellular binary and ternary IDP phase coexistence
The present theoretical development bears on the sequence dependence of multicomponent phase separation in the cell. However, because the cellular processes involve many species of biomolecules and are extremely complex [27,71], development of treatments much more elaborated than our simple theories will be needed for quantitative compar-ison with experiments. Nonetheless, it is instructive to explore whether our RPA and FH results are qualitatively consistent with what has been observed experimentally.
Of interest are fibrillarin FIB1 (323 residues) and nucleophosmin NPM1 (299 residues) from frog (Xenopus laevis) oocytes. These IDPs tend to demix, exhibiting phase behaviors that likely underpin the assembly of nucleolar subcompartments [27]. Treating histidine sidechains at pH 7 as neutral, the net charge of FIB1 is 19 and of NPM1 is −22. Their charge patterns, as quantified by SCD [Eq. (23)] = 4.126 and −0.119, respectively, are substantially different. Thus, the tendency for FIB1 and NPM1 to demix is qualitatively in line with the RPA-predicted trend in Figs. 3 and 4.
Aqueous solutions with both FIB1 and NPM1 undergo both binary and ternary liquid-liquid phase separations. In this respect, their experimental phase diagram in Fig. 4D of Feric et al. [27] is similar to our FH phase diagram in Fig. 7. The two regions of binary coexistence of one condensed (around α or around β) and one dilute (around γ) phases in Fig. 7 correspond to their "FIB1 rich/NPM1 lean" and the "NPM1 rich/FIB1 lean" areas, whereas the ternary coexistence region in Fig. 7 with two condensed (α, β) and one dilute (γ) phases [as in Fig. 1(f)] corresponds to their "3 Phase" area [27,32].
It is noteworthy that our attempts to seek numerical solutions to ternary coexistence in the RPA models studied in Fig. 3 by minimizing Eq. (28) either ended in failure or resulted in solutions with two (among three) phases essentially identical and thus reduces the solution to that of binary coexistence. Apparently, ternary coexistence requires an effective intercomponent repulsion that is substantially stronger, as in Fig. 7, than that posited by RPA. The reason is that RPA constrains the intercomponent interaction strength χ 12 to approximately the geometric mean of the two intracomponent interaction strengths χ 11 and χ 22 (Sec. 5.1) and therefore χ 12 χ 11 , χ 22 is highly unlikely if not impossible.
This consideration suggests that difference in charge pattern alone may be insufficient to account for the rather strong effective repulsion between FIB1 and NPM1, although the impact of them being not very close to being neutral remains to be investigated. (Unlike the model KE sequences in Fig. 3 with zero net charge or Ddx4 with a charge ratio = (net charge/chain length) = −1.7% [61], their charged ratios are, respectively, +5.9% and −7.4%). In addition to the electrostatic interactions among FIB1 and NPM1, other driving forces surely also contribute to their phase behaviors. For example, the presence of ribosomal RNA (rRNA) appears to be important; and the role of rRNA and inter-phase surface tensions have been modeled in lattice simulations to rationalize the FIB1/NPM1 droplet-in-droplet organization [27]. Building on our findings, much effort will be required to ascertain the precise role played by charge pattern mismatch in this intriguing phenomenon.

Cooperativity driven by concentration-dependent relative permittivity
Most analytical formulations for charged polymer solutions, including common RPA theories, treat the relative permittivity r of the solution as a constant independent of polymer concentration. However, as we noted recently [62], because of the significantly different permittivities of water ( w ≈ 80) and protein ( p ≈ 2-4) [94], the effective r of a protein solution can change dramatically with protein concentration. Indeed, proteindependent variations of dielectric properties of the aqueous medium have been shown to be relevant to globular protein stability in thermophilic species [95][96][97][98]. Because of the anticipated importance of dielectric properties to IDP energetics such as enabling a greatly enhanced propensity to phase separate [62], here we expand our theoretical exploration to two-IDP systems and consider in more detail the physical basis of various effective medium approximations that may be applied to estimate r of IDP solutions. Since the scope of this exploration is limited to establishing certain general principles, for simplicity of the discussion we let p be the relative permittivity of both IDP species in the system such that Previously we considered two models for r [62], namely the "slab model" derived by considering a planar electric capacitor, wherein Slab r which corresponds to Eq. (3) of Bragg and Pippard when the depolarizing coefficient L = 0 [99], and the Clausius-Mossotti (CM) model predicting where are proportional to the CM expression for molecular polarizability [100]. The φdependences of Slab r (φ) and CM r (φ) are very similar (Fig. 11 of [62]). Both models recognize amino acid residues and water as materials with excluded volume such that the contributions of their molecular polarizabilities to the overall dielectric property of the medium are against a vacuum background.
An alternate approach to effective medium is the Maxwell Garnett (MG) model [84] that pictures the effective dielectric property of a composite material as arising from embedding a component (IDP in our case) into an all-permeating background medium (water in our case). IDP excluded volume is not taken into account in this approach. Applying this method to our IDP solution yields where Comparing the γ p expression in Eq. (36) with Eq. (38) indicates that γ MG correspondsup to a constant-to an effective molecular polarizability of IDP material, not against vacuum but rather in a water background (γ p is mathematically equal to γ MG when w → 1, thus MG is related to CM in this respect [101]). Another approach known as the Bruggeman (BG) model [84] is an p ↔ w symmetrized form of MG. The BG φ-dependent permittivity is given by A graphical illustration of the different physical pictures assumed by the slab/CM versus the MG/BG approaches is provided by Fig. 8. Leaving aside questions about the appropriateness of their respective physical pictures for the moment, we first examine the properties of these r (φ)'s and their implications for IDP phase separation.
Variations of several properties of the r (φ)'s are shown in Fig. 9. Although all four models give r (φ = 0) = w and r (φ = 1) = p as they should, the φ-dependences of the slab/CM models are very different from that of the MG/BG models. Here we are more interested in 1/ r than r itself because 1/ r is directly proportional to Coulomb energy. Fig. 9(a) shows that 1/ r of the slab model is linear in φ, that of CM is approximately linear; but the 1/ r 's of MG and BG increase very little when φ is small and exhibit rapid increases toward the 1/ p value for φ = 1 only for φ 0.8 and 0.6, respectively. The linear and near-linear φ-dependences of the 1/ r 's for the slab and CM models are underscored by their small first and second derivatives in φ, whereas the sharp rises of 1/ r near φ ≈ 1 for the MG and BG models are illustrated by their large φ-derivatives for φ 0.7 [ Fig. 9(b) and (c)].
As described in our previous work [62], when permittivity becomes a function of IDP concentration, the last subtraction of G(k) in Eq. (3) has to be modified because G(k) is no longer linear in φ's [62]. A straightforward generalization of Eqs. (68) and (69) of Ref. [62] to the present case with two IDP sequences (but now with neither salt nor counterions) leads to the following replacement for the RPA expression in Eq. (3) to accommodate a φ-dependent r : where and, following Eq. (67) of Ref. [62], The resulting spinodal phase diagrams predicted by the four r (φ) models are shown in Fig. 10 for the (seq1, seq2) pair at T * 0 = 0.05. This temperature is chosen to facilitate comparison with constantr results because T * 0 = 0.05 corresponds to T * = 4 when w = p = 80. Fig. 10(a)-(d) shows that all four r (φ) models have large spinodal regions extending to φ 1 + φ 2 ≈ 0.8. These spinodal regions are much larger than those predicted by constantr theories [ Fig. 10(e), (f)], pointing once again to a significant cooperative effect arising from a decrease in permittivity upon IDP condensation which in turn increases electrostatic attraction and hence more IDP condensation [62]. For instance, the condensed-phase volume fractions of the slab and CM models ≈ 0.8, which represents a > 20-fold increase from the condensed-phase volume fraction of 0.03 for a constant r = w = 80 [ Fig. 10(e)]. The corresponding enhancement in the MG and BG models are even more prominent-their condensed-phase volume fraction almost reaching the φ 1 + φ 2 = 1 limit [ Fig. 10(c), (d)]-because of the sharp rises of their r (φ)'s when φ → 1 (Fig. 9). Although the precise quantitative impact of r (φ) remains to be ascertained experimentally, our theoretical results suggest strongly that φ-dependent relative permittivity should play a significant role in IDP phase separation, and that such a physical cooperative effect may help rectify some of the RPA-predicted condensedphase volume fractions [e.g. those in Fig. 3 (a)-(d)] that seem unrealistically low.
Interestingly, while the slab and CM models enlarge a single spinodal region vis-àvis that for constantr , the MG and BG models produce an additional spinodal region close to the φ-φ 2 origin. This region is similar in scope to the constantr spinodal region, and is well separated from the MG and BG models' lower boundaries at φ 1 + φ 2 ≈ 0.6 (MG) or 0.4 (BG) for their main (much larger) spinodal regions [ Fig. 10(c)-(f)]. This feature likely arises from the fact that 1/ r 's for MG and BG barely change for φ 0.2 [ Fig. 9(a)], thus the behaviors of these models at small φ's should be similar to those of constantr models. For larger φ's, however, because of the rapid increase of their 1/ r 's with φ, the MG and BG models become similar to the slab and CM models.
Although the slab/CM and MG/BG models produce similarly expanded spinodal regions for the example in Fig. 10, the difference in their predicted r (φ)'s do have important implications on the energetics of IDP phase separation. For example, the condensed-phase volume fraction of Ddx4 was recently determined to be ≈ 0.15-0.2 [86]. At φ ≈ 0.2, the slab/CM models posit a significant cooperative effect in favor of phase separation, but the MG/BG models suggest that such cooperativity is negligible unless φ 0.6 ( Fig. 9). Issues related to effective medium approximations can be intricate, as witnessed by the extensive materials-science literature on the topic [84,101]. Nevertheless, for IDP solutions, our intuition is that the slab/CM scenario is more physically plausible than the MG/BG scenario. Consistent with the slab/CM picture in Fig. 8(a), just like dissolved folded proteins, dissolved IDPs occupy volumes excluded to water. Dissolved IDPs and folded proteins have similar densities (Fig. 11) of 1.32 -1.52 g cm −3 [102] and hence similar partial molar volumes on a per-gram basis. In contrast, the MG/BG picture in Fig. 8(b) invokes a negative effective molecular polarizability for IDP that counteracts an all-permeating water medium (γ MG < 0 because p < w ). This scenario is apparently at odds with the reality of IDP excluded volume, suggesting that while the MG/BG models may be applicable in certain solid-state situations [101], they may be problematic for IDP solutions. A definitive resolution of this question awaits further theoretical and experimental investigation.

Conclusions
To recapitulate, we have taken a step to address the sequence-phase relationship in regard to how mixing/demixing of IDP components in membraneless organelles are governed by the IDPs' amino acid sequences. Going beyond mean-field FH and OV approaches, RPA provides an approximate physical account of sequence effects [61][62][63]. Our findings point to a multivalent, stochastic, "fuzzy" mode of molecular recognition in that mixing the populations of a pair of IDP sequences is favorable if their charge patterns are similar whereas population demixing is promoted when their charge patterns are very different. For the examples studied, a quantitative correlation is observed between the RPA-predicted tendency for a pair of sequences to demix in two (binary) coexisting phases and the difference in their SCD parameters. This predicted trend is qualitatively in line with the observed demixing of the nucleolar IDPs FIB1 and NPM1 because they have very different SCD's, although a comparison of the experimental FIB1/NPM1 phase diagram [27] with our RPA and FH results indicates that inclusion of non-electrostatic interactions as well as more biomolecular species in the analysis will be necessary for a quantitative theoretical account of sequence-specific ternary coexistence. It should also be noted that our current RPA formulation does not consider counterion condensation, which can be important for IDP sequences with high net charges. A recent transfer matrix theory [103] should be helpful in tackling such situations, although as it stands this theory does not address sequence dependence. Despite our theory's limitations, the simple principles of sequence dependence suggested by the present effort should already be testable by in vitro experiments on IDP polymers [70]. As illustrated by our consideration of IDP-concentration-dependent permittivity, theoretical study of IDP phase separation is only in its infancy. The logical next steps in the development of RPA theory include extending to systems with more than two sequences and sequences with non-zero net charges. Much biophysics of IDP phase separation remains to be discovered.