Charge Pattern Matching as a "Fuzzy" Mode of Molecular Recognition for the Functional Phase Separations of Intrinsically Disordered Proteins
Yi-Hsuan Lin Jacob P. Brady Julie D. Forman-Kay Hue Sun Chan
q-bio.BM



RPA theory for two charged sequences]Charge Pattern Matching 
as a “Fuzzy” Mode of Molecular Recognition for the Functional Phase
Separations of Intrinsically Disordered Proteins


Department of Biochemistry, 
University of Toronto, and Molecular
Medicine, Hospital for Sick Children, Toronto, Ontario, Canada

Departments of Molecular Genetics, Biochemistry, and Chemistry,
University of Toronto, Toronto, Ontario, Canada

Molecular Medicine, Hospital for Sick Children, 
and Department of Biochemistry, 
University of Toronto, Toronto, Ontario, Canada

Departments of Biochemistry and Molecular Genetics, 
University of Toronto, Toronto, Ontario, Canada
chan@arrhenius.med.toronto.edu




   December 30, 2023


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
(ϕ_1,ϕ_2) for the two different IDP sequences. 
The asymmetry between the two predicted coexisting phases with 
regard to their ϕ_1/ϕ_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 
ϵ_ 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 <cit.>.
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 <cit.>.
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 <cit.>
and the DEAD-box RNA helicase Ddx4 <cit.>, 
the Caenorhabditis elegans germline-specific perinuclear
RNA granules known as P-granules <cit.>
and the Ddx3 RNA helicase LAF-1 <cit.>,
as well as the stress granules triggered by integrated stress 
response <cit.> and the RNA-binding heterogeneous 
nuclear ribonucleoprotein A1 (hnRNPA1) <cit.>.
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 <cit.>, protein 
fibrillization and thus amyloid diseases <cit.>,
and various forms of neurological disorder <cit.>.



 §.§ Seeking “sequence-phase” relationships


With the advent of IDPs, the molecular biology paradigm of seeking 
“sequence-structure” relationships for globular proteins has to be
expanded to encompass “sequence-ensemble” relationships for individual 
IDPs <cit.>. 
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 <cit.>? 
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
<cit.>.
In contrast, because IDP phase separation is a multiple-chain 
property, computationally it is extremely costly to simulate using
fully atomic explicit-chain models <cit.>, notwithstanding
promising progress made by coarse-grained approaches that treat groups 
of amino acid residues of IDPs as interaction modules in continuum space
<cit.> or on lattices <cit.> and 
simulation algorithms developed recently <cit.> 
for phase separations of globular proteins <cit.>.
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 <cit.> 
is useful for describing phase separation in the biomolecular 
context <cit.>. 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 <cit.> that the Overbeek-Voorn (OV) 
theory <cit.> 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 <cit.>. The term “RPA” was first introduced by
Bohm and Pines in their quantum mechanical collective
description of electron interactions <cit.>. The analogy
of this approach with the approximate polymer theory that considers 
terms up to quadratic in particle density was recognized by
de Gennes<cit.>. As detailed elsewhere <cit.>, 
the RPA formulation we developed <cit.>, which follows largely that 
of Olvera de la Cruz<cit.>, 
is successful in providing a physical 
rationalization <cit.> 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 <cit.>. 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 <cit.>,
and that both of these properties are well correlated with the 
charge pattern parameter κ of Das and Pappu <cit.>
and the “sequence charge decoration” parameter SCD of Sawle and 
Ghosh <cit.>. In principle, these predictions are now testable 
using experimental techniques similar to those 
employed to study “IDP polymers” <cit.>.



 §.§ Sequence-dependent multiple-component IDP phase separation


Membraneless organelles are complex functional units
consisting of many protein and nucleic acid components <cit.>. 
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 <cit.>. In a similar vein, immiscible 
liquid phases have been shown to contribute to subcompartmentalization 
of the nucleolus <cit.>. 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 <cit.>. 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 <cit.>.

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 <cit.> implies that there are physical limits to 
interaction heterogeneity and structural encodability, as has been
illustrated by simple exact models <cit.>. 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” <cit.>
involving individual IDP molecules <cit.>.

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 <cit.> posit a gradual decrease in 
ϵ_ r from the ϵ_ r(ϕ=0)≈ 80 value
for pure water with increasing ϕ. But physical consideration <cit.> 
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.





§ METHODS





 §.§ 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 <cit.>.



 §.§ 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 and , 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 <cit.>.
1D ^1H 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 % ^2H_2O/90 % ^1H_2O, pH 6.5.
HEWL was purchased from BioBasic and samples 
were dissolved in 20 mM sodium citrate, 100 mM NaCl, 
10 % ^2H_2O/90 % ^1H_2O, pH 5 (lower pH was used due to 
limited solubility of HEWL in NaPi at pH 6.5).
samples were prepared according to <cit.> and dialysed 
against 20 mM NaPi, 100 mM NaCl, 5 mM tris(2-carboxyethyl)phosphine (TCEP), 
10 % ^2H_2O/90 % ^1H_2O, pH 6.5. For 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.



 §.§ Measurement of protein concentration

Protein concentrations were determined by absorbance at 280 nm (A_280) 
after dilution into 6 M guanidinium HCl, 20 mM NaPi, pH 6.5 using the 
Beer-Lambert law with extinction coefficients of 23 950, 23 950, 36 000, and 
44 309 M^-1cm^-1 for wildtype Ddx4, , HEWL, and BSA, 
respectively <cit.>.






§ 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 <cit.>. 






§ RPA THEORY FOR TWO CHARGED SEQUENCES


Based on our previous RPA formulation for a single 
sequence <cit.>, 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 <cit.>, the electric 
charges along the sequences are written as 
{σ_1}≡{σ_1^(1), σ_1^(2), σ_1^(3), 
…, σ_1^(N_1)} and 
{σ_2 }≡{σ_2^(1), σ_2^(2), σ_2^(3), 
…, σ_2^(N_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 <cit.>, 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

    f(ϕ_1,ϕ_2) ≡F_ RPAa^3/V  T 
    = -s(ϕ_1,ϕ_2)  + f_ el(ϕ_1,ϕ_2) ,

where the negative entropy -s is the entropic contribution to 
free energy in units of k_ BT. This term is given by the standard 
FH entropy of mixing for a system comprising of s_1, s_2, and solvent:

    -s(ϕ_1,ϕ_2)  = ϕ_1/N_1lnϕ_1 + 
    ϕ_2/N_2lnϕ_2 + (1-ϕ_1-ϕ_2)ln(1-ϕ_1-ϕ_2)   ,

where 1-ϕ_1-ϕ_2 is the volume fraction of solvent.
The electrostatic contribution f_ el is calculated by 
RPA <cit.>, viz.,

    f_ el(ϕ_1,ϕ_2)  = ∫_0^∞d k k^2/4π^2{ln[ 1+  G(k) ] -  G(k) },

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. <cit.> as k (i.e., kb → k), 
such that

    G(k) = 4π/k^2(1+k^2)T^*⟨ q | Ĝ_k | q ⟩,

where 4π/[k^2(1+k^2)] is from the Fourier transformation of 
Coulomb interaction with a short-range cutoff <cit.>
in units of T,

    U_ el(r) = e^2/4πϵ_0ϵ_ r T1-e^-r/b/r  ,

ϵ_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 T). Here | q ⟩ is the 
(N_1+N_2)-dimensional column vector representing the two charge 
sequences, namely q_i = σ_1^(i) for 1≤ i≤ N_1 
and q_i = σ_2^(i-N_1) for N_1+1 ≤ i≤ N_1+N_2, 
⟨ q | is the transposed row vector, 
⟨ q | Ĝ_k |q⟩≡∑_ij q_i(Ĝ_k)_ijq_j 
with (Ĝ_k)_ij 
being the i,j element of the bare two-body correlation matrix
Ĝ_k of all possible sequence-sequence correlations <cit.>,

    Ĝ_k = (
    [  Ĝ_11(k) Ĝ_12(k),;  Ĝ_21(k)  Ĝ_22(k) ]),

and Ĝ_12(k) = Ĝ_21(k). 
As in our previous studies <cit.>, we consider a simple 
formulation in which all IDPs are modeled as Gaussian chains without 
excluded volume within the RPA formalism <cit.>. In this
approximation, there is no correlation between the positions of different 
chains; hence Ĝ_12(k) = Ĝ_21(k) =0 and

    Ĝ_11(k)_ij =    ϕ_1/N_1exp( -1/6k^2|i-j|) 
    Ĝ_22(k)_ij =    ϕ_2/N_2exp( -1/6k^2|i-j|)

follow from the average of exp(i k· 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 <cit.>).
Eq. (<ref>) then becomes

    G(k) = 4π/k^2(1+k^2)T^*[ 
    ⟨σ_1 | Ĝ_11(k) | σ_1 ⟩  + 
    ⟨σ_2 | Ĝ_22(k) | σ_2 ⟩]   ,

where


    ⟨σ_1 | Ĝ_11(k) | σ_1 ⟩
    	 =    ϕ_1/N_1∑_i,j=1^N_1σ_1^(i)σ_1^(j)exp( -1/6k^2|i-j|)   ,
    
    ⟨σ_2 | Ĝ_22(k) | σ_2 ⟩
    	=    ϕ_2/N_2∑_i,j=1^N_2σ_2^(i)σ_2^(j)exp( -1/6k^2|i-j|)   .





 §.§ 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 <cit.>. The sequences are labeled here as seq1 and seq2
respectively.
The RPA free energy function for the two sequences (Fig. <ref>) 
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 
(ϕ_1^0, ϕ_2^0) 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 phase-separated state allows for a lower free energy. Accordingly,
the system undergoes phase separation to multiple coexisting phases 
with different protein volume fractions <cit.>.

The boundary between the convex and concave regions is determined by 
the saddle point condition <cit.> 

    F̂≡   
       
        = 0   .

This boundary defines a spinodal instability region in 
which F̂ < 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 
of F̂ 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 F̂ < 0
condition, viz., F̂ > 0, does not by itself
exclude the possibility that both eigenvalues of F̂ are
negative and thus the system is thermodynamically unstable despite
not satisfying F̂ < 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. (<ref>) is the 
valid condition for spinodal boundaries, examples of which are shown as 
dashed curves in Fig. <ref> for the (seq1, seq2) system.

In two-component systems such as those consisting of a single IDP species
and solvent molecules, spinodal instability necessarily leads to 
binary coexistence <cit.>. 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  <cit.> 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 <cit.>. 
As examples, three instances of binary coexistence in the 
(seq1, seq2) system are depicted in Fig. <ref> 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 <cit.>, and 
by balancing the chemical potentials for each of the components across
different phases <cit.>. 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 (ϕ_1^0, ϕ_2^0), 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 (ϕ_1^α,ϕ_2^α) and 
(ϕ_1^β,ϕ_2^β), respectively, are given by 


    f_ bulk = f(ϕ_1^0,ϕ_2^0), 
        f_ sep = v f(ϕ_1^α,ϕ_2^α) + 
    (1-v) f(ϕ_1^β,ϕ_2^β),

	

wherein conservation of volume of each of the components implies that 


    v ϕ_1^α + (1-v)ϕ_1^β    = ϕ_1^0   , 
    
    v ϕ_2^α + (1-v)ϕ_2^β    = ϕ_2^0   .

	

By rewriting Eq. (<ref>) to express IDP volume
fractions in β as functions of v with 0<v<1
and volume fractions in α:


    ϕ_1^β    = ϕ_1^0 - v ϕ_1^α/1-v, 
    
    ϕ_2^β    = ϕ_2^0 - v ϕ_2^α/1-v,


f_ sep is seen as a function
of ϕ_1^α, ϕ_2^α, and v,

    f_ sep(v,ϕ_1^α, ϕ_2^α) =
    	 v f(ϕ_1^α,ϕ_2^α) 
    	 + (1-v) f(ϕ_1^0- v ϕ_1^α/1-v,
    ϕ_2^0 - v ϕ_2^α/1-v)   .
 
Note that lim_v→ 0
f_ sep(v,ϕ_1^α, ϕ_2^α)
=lim_v→ 1f_ sep(v,ϕ_1^α, ϕ_2^α)= 
f(ϕ_1^0,ϕ_2^0), where 
(ϕ_1^β,ϕ_2^β)→ (ϕ_1^0,ϕ_2^0) for v→ 0 and
(ϕ_1^α,ϕ_2^α)→ (ϕ_1^0,ϕ_2^0) for v→ 1.

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 <cit.> using the scipy.optimize.minimize function in Scipy, a Python-based numerical package
for scientific computation <cit.>. If a given 
(v,ϕ_1^α,ϕ_2^α) is found 
to yield a minimum of f_ sep among computed f_ sep values
and also satisfies 
f_ sep(v, ϕ_1^α,ϕ_2^α) < f_ bulk(ϕ_1^0,
ϕ_2^0), 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 (ϕ_1^0, ϕ_2^0) are 
larger than f_ bulk, the bulk state (ϕ_1^0, ϕ_2^0) 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 (ϕ_1^0, ϕ_2^0) 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 <cit.>. 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

    Δ f ≡ f_ sep - f_ bulk = v f(ϕ_1^α,ϕ_2^α) 
    + (1-v) f(ϕ_1^β,ϕ_2^β) - f(ϕ_1^0,ϕ_2^0)   .

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: 

    g(ϕ_1, ϕ_2) = f(ϕ_1,ϕ_2) + a_0+a_1ϕ_1+a_2ϕ_2   .

Now, the free energy difference between phase-separated and bulk phases becomes

    Δ g ≡ g_ sep - g_ bulk = 
    Δ f + a_1 [ v ϕ_1^α + (1-v) ϕ_1^β - 
    ϕ_1^0 ] +
    a_2 [ v ϕ_2^α + (1-v)ϕ_2^β -  ϕ_2^0 ]

because a_0 in g_ sep and g_ bulk cancel. 
The two bracketed terms in Eq. (<ref>) are identically zero
because of Eq. (<ref>).
Hence Δ g = Δ f, meaning that any
linear function of ϕ's added to f, such as the
one utilized in Fig. <ref> 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 <ref>. 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. <cit.>, phase equilibrium 
implies the following equalities,


    f_1^'α    = f_1^'β
      
    f_2^'α    = f_2^'β
    
    
    μ_ w^α    = μ_ w^β
 
	

where

    f^y ≡    f(ϕ_1^y, ϕ_2^y)   , 
    
    f_x^'y≡   . 
    ∂ f(ϕ_1,ϕ_2)/∂ϕ_x|_(ϕ_1,ϕ_2)=(ϕ_1^y, ϕ_2^y), 
     
    x =     1,2    ;   y = α,β  ,

and 

    μ_ w^y ≡ f^y - ϕ_1^y f_1^'y - ϕ_2^y f_2^'y

is the chemical potential of water <cit.>. 
Making use of the volume conservation conditions in Eq. (<ref>) 
and substituting Eq. (<ref>) for ϕ_1^β and
Eq. (<ref>) for ϕ_2^β, Eq. (<ref>)
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. (<ref>).

It is straightforward to show that the set of phase-separated volume
fractions {ϕ_1^α,ϕ_2^α,ϕ_1^β,ϕ_2^β} 
determined by Eq. (<ref>) are identical to that obtained
by minimizing f_ sep in Eq. (<ref>). 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:

    J_ sep(v, ϕ_1^α, ϕ_2^α) ≡(
    ∂ f_ sep/∂ϕ_1^α
    ∂ f_ sep/∂ϕ_2^α
    ∂ f_ sep/∂ v) =  0  .
 
In other words,


    ∂ f_ sep/∂ϕ_1^α =     v ( f_1^'α - f_1^'β) = 0    ,
    	
    ∂ f_ sep/∂ϕ_2^α =     v ( f_2^'α - f_2^'β) = 0   ,
    	
    ∂ f_ sep/∂ v  =     f^α - f^β + 
    (ϕ_1^β-ϕ_1^α) f_1^'β + (ϕ_2^β-ϕ_2^α) f_2^'β = 0    ,

	
 
wherein we have utilized Eq. (<ref>)
for ϕ_1^β and Eq. (<ref>) for ϕ_2^β.
Clearly, Eqs. (<ref>) and (<ref>) 
are equivalent to Eqs. (<ref>) and (<ref>), 
respectively, and Eq (<ref>) is equivalent to 
Eq. (<ref>) by virtue of Eq. (<ref>). Q.E.D. 

Starting with (α, β) obtained by minimizing f_ sep
in Eq. (<ref>),
only those that deviate less than 0.1% from the chemical-potential-balancing 
equalities in Eq. (<ref>) 
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 <cit.>. In particular, for the set of thirty 
KE sequences of Das and Pappu <cit.> 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 <cit.> to be correlated 
with charge pattern parameters κ <cit.> and

    SCD≡1/N∑_i=1^N∑_j=i+1^N σ_i σ_j 
    √(j-i)  ,

where i,j label the residues with charges σ_i,σ_j along a chain 
of length N <cit.>. The κ and SCD parameters exhibit 
similar correlations with single-chain radius of gyration 
R_ g <cit.>. 
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 <cit.>.
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. (<ref>), 

    ln[1 +  G(k)] -  G(k) =     -1/2 G(k)^2 + 1/3 G(k)^3 + ... 
    
    =     -1/2⟨σ_1 | Ĝ_11^*(k) | σ_1 ⟩^2 
    -1/2⟨σ_2 | Ĝ_22^*(k) | σ_2 ⟩^2  
        - ⟨σ_1 | Ĝ_11^*(k) | σ_1 ⟩⟨σ_2 | Ĝ_22^*(k) | σ_2 ⟩
        + O( G(k)^3)   ,

where Ĝ_11^*(k) and Ĝ_22^*(k)
are the product of 4π/[k^2(1+k^2)T^*] with, respectively, the
Ĝ_11(k) and Ĝ_22(k) in Eq. (<ref>).
The first two terms after the second equality in Eq. (<ref>) are
self-interactions of the two sequences, identical to those in one-sequence 
RPA theory (see, e.g. Eq. (1) in Ref. <cit.>). 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 <ref> 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. <ref>
are computed at the same reduced temperature T^*=4. Noting that
condensed-phase 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 <ref>. T^*=4 is much higher than 
the T^* = 0.55 equivalent of room temperature (T = 300 K) when an aqueous
ϵ_ r=80 is assumed <cit.>. 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. <ref>, the area 
surrounded by blue dots and “shaded” by black dashed lines
is the region of binary coexistence. In other words,
bulk-state 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 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. <ref>(a), (b)] to negative for 
very differently patterned sequences [Fig. <ref>(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 ϕ_1^0→ 0 or ϕ_2^0→ 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. <ref>(f). As shown
in Fig. <ref>(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. <ref>(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. <ref>(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^β≃ϕ_2^0 [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. We first
consider a rather intuitive measure 

    Δ_αβ( ϕ_1/ϕ_2 ) ≡⟨| ϕ_1^α/ϕ_2^α - ϕ_1^β/ϕ_2^β| 
    ⟩

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. <ref>.
Therefore, we also consider another composition asymmetry measure 

    A_αβ≡⟨2/π| 
    tan^-1(ϕ_1^α/ϕ_2^α) - 
    tan^-1(ϕ_1^β/ϕ_2^β) | ⟩

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. <ref> shows the
variation of Δ_αβ(ϕ_1/ϕ_2) as well as
A_αβ with the difference in SCD values
of the six sequence pairs in Fig. <ref>. 
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 <cit.>. 
Accordingly, FH χ interaction parameters
for IDP sequences have to be provided phenomenologically by experiment 
or theoretically by microscopic physical 
theory.[In the caption describing the FH results in Fig. 8
of Ref. <cit.>, r_d is in fact the symbol r for
the equilibrium spacing in Eq.(S19) of Ref. <cit.>. This
typographical error does not affect the results. 
It should also be noted that because Ref. <cit.> 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. <cit.>, however,
the effective Debye length in Fig. 8 of Ref. <cit.> was also 
set to 4.3 Å.] 
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. (<ref>).
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. (<ref>)
because G(k)∝ϕ [Eqs. (<ref>) and (<ref>)]. 
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

    f_ int^ FH = -( χ_11ϕ_1^2 + χ_22ϕ_2^2 + 
    2χ_12ϕ_1ϕ_2 )

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. <ref>.
If we identify the ∫ dk k^2/4π^2 integral [Eq. (<ref>)] 
of the order G(k)^2 terms in the RPA expansion
Eq. (<ref>) with the FH interaction term in 
Eq. (<ref>), 
we may define
χ̃_11(k)ϕ_1^2≡
2{⟨σ_1 | Ĝ_11(k)|σ_1⟩/[ k (1+k^2) T^*]}^2
and
χ̃_22(k)ϕ_2^2≡
2{⟨σ_2 | Ĝ_22(k)|σ_2⟩/[ k (1+k^2)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).

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. <ref>) exhibit a trend similar
to that seen in the RPA-predicted Fig. <ref>.
Specifically, Fig. <ref>(a) is similar to 
Fig. <ref>(a), Fig. <ref>(b) 
to Fig. <ref>(c), and Fig. <ref>(c) 
to Fig. <ref>(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. <ref>(a) 
Fig. <ref>(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. <ref>(c) with that in
Fig. <ref>(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. (<ref>) 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. <ref> 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. <ref>(a), it is clear that
a stronger inter-component attraction (larger χ_12)
makes the phase-separated region bulge [Fig. <ref>(a)],
whereas a weaker inter-component attraction (smaller χ_12)
shrinks it [Fig. <ref>(b)]. 
Nonetheless, the tie-line slopes are positive in both situations, indicating
that the two components are largely miscible.

Fig. <ref>(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 ϕ_1 ≪ϕ_2
[Fig. <ref>(c), shaded area and inset].
The tie lines in these two regions are almost parallel to either
the ϕ_1- or the ϕ_2-axis, implying that one
component is dominant while the concentration of the other component
barely changes upon phase separation. 
With such an effective inter-component repulsion
(i.e., less favorable inter-component attraction
vis-à-vis the strengths of intra-component cohesion), an 
additional phase-separated regime of poor miscibility similar
to that in Figs. <ref>(f) and <ref>(c)
is induced, wherein all tie-line slopes are negative, 
signaling substantial demixing
[Fig. <ref>(c), region close to ϕ_1+ϕ_2=1 
with tie-line slopes =-1]. 



 §.§ Ternary coexistence in FH model


If χ_12 is made even weaker than that in 
Fig. <ref>(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 <cit.>. 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. <ref>(c) as
starting point. Here the three phase-separated regions are close to
the three boundaries, and their tie lines are essentially parallel to
the respective boundaries. Under this circumstance, when effective
inter-component 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. <ref>). 

In order to determine the three phases in ternary coexistence mathematically, 
we extend the phase-separated free energy expression
in Eq. (<ref>) for phases (α,β) to including one 
additional phase γ, viz.,

    f_ ternary = 
    v_α f(ϕ_1^α,ϕ_2^α) 
    + v_β f(ϕ_1^β,ϕ_2^β)
    + (1-v_α-v_β) f(ϕ_1^γ,ϕ_2^γ)   ,

where the fractional volumes v_α, v_β are functions of 
the six ϕ's by virtue of volume conservation,

    ∑_y v_y ϕ_x^y + (1- ∑_y v_y) ϕ_x^γ=ϕ_x^0,

where x=1,2 and ∑_y is over y=α,β. 
Now the equalities in Eq. (<ref>) have to include the 
addition phase to become


    f_1^'α    = f_1^'β = f_1^'γ  ,  
      
    f_2^'α    = f_2^'β = f_2^'γ  ,   
    μ_ w^α    = μ_ w^β = μ_ w^γ  .
 
	

Because here we have six equalities in Eq. (<ref>) 
for six phase-separated ϕ_x^y'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 (ϕ_1^0,ϕ_2^0)'s when they are
on different tie lines.

Similar to the binary coexistence case in Sec. <ref>,
we proceed to demonstrate that minimizing Eq. (<ref>) 
is equivalent to solving the equations in Eq. (<ref>). 
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. (<ref>) to express ϕ_1^γ 
and ϕ_2^γ as

    ϕ_x^γ = ϕ_x^0-_y v_y ϕ_x^y/1- _y v_y  ,

where x,y and ∑_y have the same meanings as above. In this
notation, Eq. (<ref>) becomes

    f_ ternary 
    (v_α,v_β,ϕ_1^α,ϕ_2^α,ϕ_1^β,ϕ_2^β)
    = ∑_y v_y f^y + (1-∑_y v_y) f^γ  .

Similarly to Eq. (<ref>) for binary phase separation, we 
calculate the six derivatives of f_ ternary and set them
to zero as in Eq. (<ref>) as necessary conditions for
the minimization of f_ ternary, resulting in


    ∂ f_ ternary/∂ϕ_x^y =     v_y f^' y_x 
    - (1-∑_y^' v_y^') 
    f_x^'γ·-v_y/ 1-_y^' v_y^' 
    = v_y ( f^' y_x - f_x^'γ) = 0    ,
    
    ∂ f_ ternary/∂ v_y =     f^y - f^γ 
    + (1-∑_y^' v_y^')  
    [ ∑_x f_x^'γ·-ϕ_x^y(1-_y^' v_y^') +ϕ_y^0 
    - _y^' v_y^'ϕ_x^y^'/(1-_y^' v_y^') ^2] 
    
    
    =     f^y -f^γ + ∑_x f^' γ_x (ϕ^γ_x - ϕ^y_x) 
    
    
    =    μ_w^y - μ_w^γ = 0   ,


where ∑_y^' sums over y^'=α,β.
Substituting x=1,2 in Eq. (<ref>) yields 
Eqs. (<ref>a) and (<ref>b), 
whereas substituting y=α, β in Eq. (<ref>) yields 
Eq. (<ref>c). Q.E.D.





Fig. <ref> 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. <ref>(a) and turquoise lines 
in Fig. <ref>(b) and (c)]. The thermodynamic stability
of the ternary phase-separated state within this region is illustrated
by the Δ f quantity plotted in Fig. <ref>(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. <ref>(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 (ϕ_1^0,ϕ_2^0) in the ternary
region, fractional volumes v_α, v_β, 
and v_γ in the respective coexisting phases α, β, 
and γ are determined by solving
Eq. <ref> and setting v_γ=1-v_α-v_β.
In terms of the three-dimensional vectors
Φ_0≡(ϕ_1^0,ϕ_2^0,0),
Φ_z≡(ϕ_1^z,ϕ_2^z,0),
Φ_0z≡Φ_z-Φ_0,
and Φ_z_1z_2≡Φ_z_2-Φ_z_1
where z,z_1,z_2=α,β,γ, 
v_α=|Φ_0γ×Φ_βγ|
/|Φ_αγ×Φ_βγ|,
v_β=|Φ_0α×Φ_αγ|
/|Φ_αγ×Φ_βγ|
and
v_γ=|Φ_0β×Φ_αβ|
/|Φ_αγ×Φ_βγ|.
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. <ref>.




§ 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 <cit.>, development of treatments much
more elaborated than our simple theories will be needed for quantitative 
comparison 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 <cit.>. 
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. (<ref>)] =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. <ref> and <ref>.

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. <cit.> is similar
to our FH phase diagram in Fig. <ref>. The
two regions of binary coexistence of one condensed (around α or
around β) and one dilute (around γ) phases in 
Fig. <ref> correspond to their “FIB1 rich/NPM1 lean” 
and the “NPM1 rich/FIB1 lean” areas,
whereas the ternary coexistence region in Fig. <ref>
with two condensed (α,β) and one dilute (γ) phases
[as in Fig. 1(f)] corresponds to their “3 Phase” 
area <cit.>.

It is noteworthy that our attempts to seek numerical solutions 
to ternary coexistence in the RPA models studied in Fig. <ref> 
by minimizing Eq. (<ref>) 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. <ref>, 
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. <ref>) 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. <ref> with zero net charge
or Ddx4 with a charge ratio = (net charge/chain length) 
=-1.7% <cit.>,
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 <cit.>.
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 <cit.>, because of the significantly different 
permittivities of water (≈ 80) and protein 
(≈ 2–4) <cit.>, the effective ϵ_ r of
a protein solution can change dramatically with protein concentration. 
Indeed, protein-dependent variations of dielectric properties of the 
aqueous medium have been shown to be relevant to globular protein stability
in thermophilic species <cit.>. 
Because of the anticipated importance of dielectric properties to IDP 
energetics such as enabling a greatly enhanced propensity to phase 
separate <cit.>, 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 
be the relative permittivity of both IDP species in the system such that 
ϵ_ r(ϕ_1,ϕ_2)=ϵ_ r(ϕ) where
ϕ=ϕ_1+ϕ_2.

Previously we considered two models for ϵ_ r <cit.>,
namely the “slab model” derived by considering a planar electric
capacitor, wherein

    ϵ_ r^ Slab(ϕ) = /ϕ + (1-ϕ)  ,

which corresponds to Eq. (3) of Bragg and Pippard when the depolarizing 
coefficient L=0 <cit.>, and the Clausius-Mossotti (CM) model
predicting

    ϵ_ r^ CM(ϕ) = 1+2[(1-ϕ) + 
    ϕ ]/1-[(1-ϕ) + ϕ ]  ,

where

    ≡ -1/ + 2  ,
        ≡ -1/ + 2

are proportional to the CM expression for molecular 
polarizability <cit.>. The ϕ-dependences of 
ϵ_ r^ Slab(ϕ) and ϵ_ r^ CM(ϕ)
are very similar (Fig. 11 of <cit.>). 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 <cit.> 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

    ϵ_r^ MG(ϕ) = 
    1+2ϕγ_ MG/1-ϕγ_ MG,

where

    γ_ MG≡ -/ +2  .

Comparing the γ_ p expression in Eq. (<ref>) with 
Eq. (<ref>) indicates
that γ_ MG corresponds—up 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 → 1,
thus MG is related to CM in this respect <cit.>).
Another approach known as the Bruggeman (BG) model <cit.>
is an ↔ symmetrized form of MG. The BG ϕ-dependent
permittivity is given by

    ϵ_r^ BG(ϕ) ≡b_ BG(ϕ) + 
    √(8+b_ BG(ϕ)^2)/4  ,

where

    b_ BG(ϕ) = [2ϕ-(1-ϕ)] + [2(1-ϕ)-ϕ]  .

A graphical illustration of the different physical pictures assumed by 
the slab/CM versus the MG/BG approaches is provided by
Fig. <ref>.
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. <ref>.  Although all four models 
give ϵ_ r(ϕ=0)= and ϵ_ r(ϕ=1)= 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. <ref>(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/
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. <ref>(b) and (c)].




As described in our previous work <cit.>, when permittivity 
becomes a function of IDP concentration, the last subtraction of 
G(k) in Eq. (<ref>) has to be modified
because G(k) is no longer linear in ϕ's <cit.>. 
A straightforward generalization of Eqs. (68) and (69) of Ref. <cit.>
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. (<ref>) to accommodate a ϕ-dependent 
ϵ_ r:

    f_ el(ϕ_1,ϕ_2)  = ∫_0^∞d k k^2/4π^2{ln[ 1+  G_1(k) ] -  G_2 (k) }  ,

where

    G_1(k) =    4π/k^2(1+k^2)T^*_0ϵ_r(ϕ)[ 
    ⟨σ_1 | Ĝ_11(k) | σ_1 ⟩  + 
    ⟨σ_2 | Ĝ_22(k) | σ_2 ⟩], 
     G_2(k) =    4π/k^2(1+k^2)T^*_0ϵ_r(ϕ)[ 
    ϕ_1/N_1∑_i| σ_1^(i) | 
    	+ ϕ_2/N_2∑_i| σ_2^(i) | ]   ,

and, following Eq. (67) of Ref. <cit.>,

    T_0^* ≡ 4πϵ_0 k_ BTb/e^2 = T^*/ϵ_ r  .

The resulting spinodal phase diagrams predicted by the four 
ϵ_ r(ϕ) models are shown in 
Fig. <ref> for 
the (seq1, seq2) pair at T_0^* = 0.05. This temperature is chosen
to facilitate comparison with constant-ϵ_ r results
because T_0^* = 0.05 corresponds to T^*= 4 when ==80.  






Fig. <ref>(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 constant-ϵ_ r theories
[Fig. <ref>(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 <cit.>.
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= =80 [Fig. <ref>(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. <ref>(c),
(d)]—because of the sharp rises of their ϵ_ r(ϕ)'s 
when ϕ→ 1 (Fig. <ref>). 
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 condensed-phase volume fractions 
[e.g. those in Fig. <ref> (a)–(d)] that seem 
unrealistically low.

Interestingly, while the slab and CM models enlarge a single spinodal region
vis-à-vis that for constant-ϵ_ r,
the MG and BG models produce an additional spinodal region
close to the ϕ-ϕ_2 origin. This region is similar in scope to 
the constant-ϵ_ r 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. <ref>(c)–(f)].
This feature likely arises from the fact that 
1/ϵ_ r's for MG and BG barely change
for ϕ≲ 0.2 [Fig. <ref>(a)],
thus the behaviors of these models at small ϕ's should be similar
to those of constant-ϵ_ r 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. <ref>,
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 <cit.>.
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. <ref>). Issues related to effective
medium approximations can be intricate, as witnessed by the extensive 
materials-science literature on the topic <cit.>.
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. <ref>(a), just like
dissolved folded proteins, dissolved IDPs occupy volumes excluded to water. 
Dissolved IDPs and folded proteins have similar densities 
(Fig. <ref>) of
1.32 – 1.52 g cm^-3 <cit.> and hence
similar partial molar volumes on a per-gram basis.
In contrast, the MG/BG picture in Fig. <ref>(b) 
invokes a negative effective molecular polarizability 
for IDP that counteracts an all-permeating water medium 
(γ_ MG<0 because <). 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 <cit.>, 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 <cit.>.
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 <cit.> 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 <cit.> 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 <cit.>.
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.






§ ACKNOWLEDGEMENTS

We thank Alaji Bah, Robert Vernon (Hospital for Sick Children), 
Lewis Kay (University of Toronto), Che-Ting Chan (Hong Kong University 
of Science and Technology), Kingshuk Ghosh (University of Denver), 
Pak-Ming Hui (Chinese University of Hong Kong), and Huan-Xiang Zhou
(University of Illinois at Chicago) for helpful 
discussions. We are also grateful to an anonymous referee for
suggestions that led to the pedagogical graphics in Fig. 1. 
This work was supported by Canadian Cancer Society Research 
Institute grant no. 703477, Canadian Institutes of Health Research grant 
MOP-84281, and computational resources provided by SciNet of 
Compute/Calcul Canada.







§ REFERENCES

 
iopnum


