August 29, 2022

Analytical Formulation and Field-Theoretic

Simulation of Sequence-Specific Phase Separation of

Proteinlike Heteropolymers with Short- and

Long-Spatial-Range Interactions

Jonas WESSÉN, Suman DAS, Tanmoy PAL, and Hue Sun CHAN∗

Department of Biochemistry, University of Toronto, Toronto, Ontario M5S 1A8,
Canada

$*$Corresponding author
† E-mail<EMAIL_ADDRESS>Tel: (416)978-2697; Fax: (416)978-8548
† URL: http://biochemistry.utoronto.ca/person/hue-sun-chan/
† Mailing address:
† Department of Biochemistry, University of Toronto, Medical Sciences Building
-- 5th Fl.,
† 1 King’s College Circle, Toronto, Ontario M5S 1A8, Canada.

Abstract

A theory for sequence dependent liquid-liquid phase separation (LLPS) of
intrinsically disordered proteins (IDPs) in the study of biomolecular
condensates is formulated by extending the random phase approximation (RPA)
and field-theoretic simulation (FTS) of heteropolymers with spatially long-
range Coulomb interactions to include the fundamental effects of short-range,
hydrophobic-like interactions between amino acid residues. To this end, short-
range effects are modeled by Yukawa interactions between multiple
nonelectrostatic charges derived from an eigenvalue decomposition of pairwise
residue-residue contact energies. Chain excluded volume is afforded by
incompressibility constraints. A mean-field approximation leads to an
effective Flory $\chi$ parameter, which, in conjunction with RPA, accounts for
the contact-interaction effects of amino acid composition and the full
sequence effects of long-range electrostatics in IDP LLPS, whereas FTS based
on the formulation provides full sequence dependence for both short- and long-
range interactions. This general approach is illustrated here by applications
to variants of a natural IDP in the context of several different amino-acid
interaction schemes as well as model hydrophobic-polar sequences sharing the
same composition. Effectiveness of the methodology is further verified by
coarse-grained explicit-chain molecular dynamics simulations.

INTRODUCTION

Tremendous recent advances have made it abundantly clear that biomolecular
condensates serve many important biological functions.rosen2017 ; cliff2017 ;
rosen2021 The architecture of functional intraorganismic biomolecular
condensates can be highly complex. Some contain hundreds of different kinds of
proteins---including folded protein domains, intrinsically disordered proteins
(IDPs) and intrinsically disordered protein regions (IDRs)---and various
nucleic acids participating in an intricate network of interactions, as
exemplified by the $\sim 300$ protein species identified from proteomic
analyses of yeast and mammalian stress granules using mass
spectrometry.parker2016 While the existence of intracellular liquid-like
compartments has been suggested more than 120 years ago from observing that
the protoplasm of starfish and sea urchin eggs resembles ‘‘a mixture of
liquids’’ with suspended droplets of different chemical nature,wilson1899 the
intense modern interest in biomolecular condensates was sparked by the
recognition that the assembly of these condensates bears close resemblance to
the physical process of liquid-liquid phase separation (LLPS)---as noted in a
seminal study of germline P granules a little more than a decade
ago.brangwynne2009 The LLPS perspective has since contributed greatly to the
discovery and understanding of a large variety of biomolecular condensates.
Progress is often achieved through experimental characterizations of the phase
behaviors as well as the biochemical and materials properties of constructs,
composing of only a few species of protein and nucleic acid, to serve as
models of the biomolecular condensates of interest.Rosen12 ; McKnight12 ;
Nott15 ; tanja2015 ; cliff2015 ; parker2015 ; Michnick2016 ; babu2018 ;
Monika2018Rev ; Roland2019 ; shorter2019 One novel observation from these
investigations is that multivalent, relatively weak interactions of IDPs and
IDRs often play a major role in the formation of biomolecular
condensates.McKnight12 ; Nott15 ; tanja2015 ; parker2015 These interactions
among largely disordered conformations of IDPs and IDRs are transient and
presumably nonstoichiometric,Fawzi2015 ; jacob2017 though in some instances
they entail sampling disordered chain configurations with labile but
nonetheless specific fibril-like local structures McKnight12 ; McKnight2018 .
Despite the stochastic nature of these IDP/IDR interactions, they are capable
of contributing to a ‘‘fuzzy’’, sequence-dependent mechanism of molecular
recognition.feric2016 ; njp2017

Confronted by the overwhelming complexity of functional intraorganismic
condensates, it goes without saying that few-component LLPS under
thermodynamic equilibrium is a rudimentary---yet important---model for gaining
insights into biomolecular condensates in living organisms. For instance,
nonequilibrium behaviors of in vivo condensates are more appropriately viewed
as those of active liquids;julicher2018 ; lee2022 and formulations beyond
simple LLPS theory are needed to describe the size, viscoelastic, and other
materials properties of biomolecular condensates.dufresne2020 ; hxzhou2021 ;
espinosa2021 In addition to physio-chemical processes akin to LLPS, gelation
and percolationRohitRosen2017 ; biochemrev as well as structure-specific
stoichiometric interactionsTjian2019 ; musacchio2022 are expected to also
contribute significantly to the formation of biomolecular condensates. It is
likely that these processes are coupledRohit2022 as, at least in the example
of a model condensate for postsynaptic densities, stoichiometric interactions
alone do not account for the experimentally observed phase
properties.YHLin2022 Nonetheless, despite these complexities, since LLPS is
integral to---though not exclusively responsible for---biomolecular condensate
formation, experimentalbrangwynne2009 ; Rosen12 ; McKnight12 ; Nott15 ;
tanja2015 ; cliff2015 ; parker2015 ; Michnick2016 ; babu2018 ; Monika2018Rev ;
Roland2019 ; shorter2019 and theoretical/computationalnjp2017 ; CellBiol ;
NatPhys ; linPRL ; linJML ; lin2017 ; dignon18 ; suman1 ; jeetainPNAS ; suman2
; stefan2019 ; joanElife ; joanJPCL studies of equilibrium LLPS is a valuable
tool for gaining physical, chemical, and biological insights. In this regard,
they play an essential role similar to other simple model systems in many
branches of science---including, as a befitting example for this Special
Issue, a seminal lattice model study by Onuchic and coworkersleopold that
inspired the funnel landscape picture of protein folding,wolynes ; dillchan97
with the obvious proviso that many aspects of the recent simple model LLPS
systems’ quantitative relationship with the biological functions of
intraorganismic biomolecular condensates remain to be delineated. Indeed, the
LLPS perspective has already led to several notable conceptual advances,
including a likely role of the general physical principle of phase separation
in biomolecular homeostasis,njp2017 ; biochemrev ; Zechner2020 ; safran2021
novel clues to neurological effects of hydrostatic pressure,roland2020 ;
roland2022 and a likely link between neurological disease-causing mutations
and bioinformatics-inferred LLPS propensities,BrianJulie2020 to name a few.

Depending on the phenomena and scientific questions of interest,
theoretical/computational modeling of biomolecular LLPS may endeavor to
capture different levels of structural and energetic detail.biochemrev These
modeling approaches include---but certainly not limited to---mean-field Flory-
Huggins (FH)FH-ref and Overbeek-Voorn (OV)OV-ref theories of polymer
solutions, random phase approximation (RPA) theories to address effects of
sequence charge patterns on LLPS of polyampholytes as model IDPs/IDRs,linPRL ;
linJML ; delacruz2003 field-theoretic simulation (FTS)Fredrickson2006 to
provide an improved account of charge and matter density fluctuations beyond
RPA for polyampholyte LLPSjoanElife ; joanJPCL ; joanPNAS ; joanJCP \---though
FTS is still limited by finite-size effects and in its treatment of excluded
volumeirback2020 ; Pal2021 ; irback2021 ; irback2022 \---as well as explicit-
chain simulations of lattice,feric2016 ; suman1 ; stefan2019 ; lassi2019
continuum coarse-grained,dignon18 ; jeetainPNAS ; suman2 ; SumanPNAS ;
panag2017 ; jeetain-rev2021 and atomicregis2017 ; JeetainAtom models of
IDP/IDR LLPS. Transfer matrix, restricted primitive model simulations, and
other techniques have been applied to study complex coacervation of
polyampholytes and polyelectrolytes;transferM2017 ; singperry2017 ; sing2020
and LLPS of single or multiple biomolecular species, especially those
involving folded proteins and folded domains, have also been modeled by patchy
particlesvlachy2016 ; hxzhou2018 ; roxana2020 using simulationshxzhou2018 ;
roxana2020 as well as analytical formalisms based on Wertheim’s thermodynamic
perturbation theory.vlachy2016 ; wertheim1986 As expected, the required
computation increases with structural and energetic details a model seeks to
capture. In this context, it is noteworthy that even basic FH theory, which
requires minimal numerical effortmimb2022 that has recently been made even
more tractable by a novel self-consistent solution,analyticalFH can be very
useful in advancing knowledge about biomolecular LLPS. This is exemplified by
the applications of FH to delineate scenarios of tie-line patterns and their
ramifications for homeostasis,njp2017 ; biochemrev ; safran2021 to ascertain
the extent of void-volume contributions to the hydrostatic pressure dependence
of LLPS,roland2020 and to rationalize experimental data on the impact of
aromatic valence on IDP LLPS.TanjaScience2020

A fundamental limitation of mean-field FH and OV theories is that they
consider the composition of heteropolymer sequences without accounting for the
full effects of the sequential arrangements of monomers (e.g., amino acid
residues for IDPs/IDRs) along heteropolymer chains. While explicit-chain
simulations of LLPS readily embody sequence-pattern effects,dignon18 ; suman1
; jeetainPNAS ; suman2 ; stefan2019 ; regis2017 ; SumanPNAS ; JeetainAtom ;
Alan2020 ; koby2020 ; koby2022 theories that take into account sequence
patterns but are computationally less intensive have proven useful as
complementary approaches, especially for screening large number of sequences.
To date, these sequence-specific theories, which include RPA, FTS, variational
approaches,kings2015 and RPA augmented by Kuhn-length renormalizationrG-RPA-
ref or a more physical account of charge density and its
fluctuations,dePablo2022 have focused on spatially long-range Coulomb
interactions entailed by the heteropolymeric sequence patterns of electric
chargeskings2015 ; rohit2013 ; Kings2022 and their effects on LLPSlinPRL ;
linJML ; lin2017 ; joanElife ; joanJPCL ; joanPNAS ; joanJCP ; Pal2021 ;
mimb2022 ; rG-RPA-ref ; wessen2021 ; wessen2022 as well as single-chain
polyampholyte/IDP conformational propertieslin2017 ; dePablo2022 ; Kings2022 ;
kings2017 ; firman2018 ; huihui2018 ; kings2020 and the interaction of a pair
of polyampholyte chains.Alan2020

These theoretical treatments of spatially long-range electrostatic
interactions have led to many important new physical insights; but these
theories by themselves do not address the full sequence-dependent effects of
spatially short-range, contact-like hydrophobicregis2017 and
$\pi$-relatedSongKAW2013 ; robert interactions, which are integral parts of
the physical driving forces in biomolecular LLPSdignon18 ; SumanPNAS ;
regis2017 ; robert ; moleculargrammar ; kitahara2021 ; kameda2022 and IDP
conformational propertiesjeetainPNAS ; zhengHP ; song21 and
interactions.SongKAW2013 In lieu of an account of sequence-pattern dependence,
short-range interactions are sometimes treated by augmenting a sequence-
dependent theory for electrostatic interactions such as RPA by a composition-
dependent mean-field FH account of spatially short-range interactions.linPRL ;
linJML One source of difficulty in extending RPA to include spatially short-
range attractive (favorable) interactions is that they are prone to produce
mathematical singularities in the RPA formulation. This observation led us to
consider the possibility of incorporating these interactions in FTS because
FTS does not require the approximations in RPA theory and should, therefore,
avoid the unphysical singularities. We put forth such a theory here. To
account for nonelectrostatic, spatially short-range interactions in our FTS
formulation, additional fields are introduced, wherein different amino acid
residues are assigned different nonelectrostatic charges associated with the
newly introduced fields. The nonelectrostatic charges are determined as
components of eigenvectors using eigenvalue decomposition of the pairwise
interaction energiesHaoLi1997 ; chan99 from either the values for the 20 amino
acids on a hydrophobicity (hydropathy) scaledignon18 ; Rossky ; Urry ; Urry-
Mittal ; FB or the $20\times 21/2=210$ pairwise contact energies between the
amino acids (the interaction matrix) derived from statistical analyses of the
Protein Data Bank (PDB) and/or from physical considerations.dignon18 ; MJ85 ;
MJ96 ; KH ; Mpipi ; EDLevy2022 This approach is computationally tractable
because matrices for physical pairwise interactions among amino acid residues
tend to possess a dominant eigenvector.HaoLi1997 ; chan99 ; godzik1995 ;
chanMIT Thus, only a small number of types of nonelectric charge ($\ll 20$)
are necessary to provide a mathematically good approximation to the $210$
pairwise interactions. In view of the long-standing recognition that there are
significant discrepancies among hydrophobicity scales because of differences
in the experimental or computational techniques used for their
construction,pakarplus1997 ; devido1998 ; Pappu2021 the size- and length-scale
dependence of hydrophobic interactionschandler1999 ; scheraga2007 ;
scheraga2008 and their deviations from pairwise additivitychan2011 ; shimizu01
; shimizu02 (as manifested, e.g., by the fact that pairwise interactions alone
are insufficient to account for protein folding cooperativitychan2011 ;
Chan2004 ; chan98 ), we apply our new formulation to several different
hydrophobicity scalesdignon18 ; Urry-Mittal ; FB and pairwise interaction
schemesdignon18 ; Mpipi to compare their predicted LLPS behaviors for the
wildtype and variant IDRs of the DEAD-box RNA helicase Ddx4jacob2017 ; Nott15
and model two-letter (hydrophobic-polar) sequences.Statt2020 ; panagio2021 The
theoretical predictions are compared against experimental trend and
corresponding predictions from coarse-grained molecular simulations. These
results together with their underpinning theoretical development and their
ramifications for future efforts are provided in detail below.

MODELS AND METHODS

We begin by presenting a general analytical formulation, amenable to RPA and
FTS, for a model system containing multiple chains of heteropolymers with
specific sequences of residues (monomers) together with small ions and solvent
molecules, wherein the constituents may interact via spatially short-range,
contact-like interactions as well as spatially long-range Coulomb
interactions.

Spectral (Eigenvalue) Decomposition of Residue-Residue Energy Matrices and
Definition of Nonelectric Charges. Assuming that the spatially short-range
interactions among the twenty types of amino acid residues may be described
approximately by a $20\times 20$ matrix $\varepsilon_{r,r^{\prime}}$ of
contact energies for the pairwise interactions between residue types $r$ and
$r^{\prime}$, this matrix can be written in the diagonalized form

$\varepsilon_{r,r^{\prime}}=\sum_{a=1}^{20}\lambda_{a}q^{(a)}{}_{r}q^{(a)}{}_{r^{\prime}}$
(1)

because $\varepsilon_{r,r^{\prime}}$ is symmetric
($\varepsilon_{r,r^{\prime}}=\varepsilon_{r^{\prime},r}$), with
$\\{\lambda_{a}\\}_{a=1}^{20}$ being the eigenvalues of
$\varepsilon_{r,r^{\prime}}$, and $q^{(a)}{}_{r}$ constituting a set of twenty
orthogonal, unit-normalized eigenvectors labeled by $a$. We may then view the
eigenvector component $q^{(a)}{}_{r}$ as the value of a mathematically defined
type-$a$ nonelectric charge for residue type $r$. Without loss of generality,
we label the eigenvalues (and thus the nonelectric charge type) such that the
eigenvalues are ordered according to descending magnitude, i.e.,
$|\lambda_{a}|\geq|\lambda_{a+1}|$.

The spectral form in Eq. (1) implies that we can decompose the full
$\varepsilon_{r,r^{\prime}}$ into a set of interactions, with overall
strengths $\lambda_{a}$ and relative residue-dependent strengths given by
product of nonelectric charges $q^{(a)}{}_{r}q^{(a)}{}_{r^{\prime}}$. Indeed,
this approach has long been used to construct hydrophobicity scales using the
eigenvalue decomposition of statistical contact potentials such as the classic
Miyazawa-Jernigan (MJ) matrix.HaoLi1997 ; chan99 ; MJ85 ; MJ96

Although all twenty eigenvalues are needed in Eq. (1) to recover the full
$\varepsilon_{r,r^{\prime}}$ exactly, one can often achieve a very good
numerical agreement by including only the first few most dominant
$\lambda_{a}$, i.e. by truncating the sum at some $\tilde{a}\leq 20$,

$\sum_{a=1}^{20}\rightarrow\sum_{a=1}^{\tilde{a}},$ (2)

because of the pattern of hydrophobic-like interactions among amino acid
residues.HaoLi1997 ; chan99 ; godzik1995 ; chanMIT For our present purpose, we
consider two interaction matrices proposed recently for the study of IDP LLPS,
namely, $\varepsilon^{\rm KH}_{r,r^{\prime}}$ obtained through a shifted MJ
matrix,dignon18 ; MJ96 and $\varepsilon^{\rm Mpipi}_{r,r^{\prime}}$, which is
partly based on explicit-water simulations of pairwise amino acid residue
interactions and has recently been used to successfully predict LLPS
propensities of several biologically important IDPsMpipi (Fig. 1). The
accuracies of truncated eigenvalue decompositions of $\varepsilon^{\rm
KH}_{r,r^{\prime}}$ and $\varepsilon^{\rm Mpipi}_{r,r^{\prime}}$ [Eq. (2),
with $\tilde{a}=1,2,3,4,$ and $5$] are illustrated in Fig. 1a,b and the
associated nonelectric charges associated with the $a=1,2,3,4$, and $5$
eigenvalues are provided in Fig. 1c.

Fig. 1: Eigenvalue decomposition of (a) $\varepsilon^{\rm KH}_{r,r^{\prime}}$
and (b) $\varepsilon^{\rm Mpipi}_{r,r^{\prime}}$ using Eqs. (1) and (2) where
the summations over eigenvalues are truncated at $\tilde{a}$. Accuracy of a
set of truncated summations (vertical variable) is indicated by the proximity
of the scatter plot [red (a) and green (b) data points] to the slope $=1$
black line through the origin. (c) Nonelectric charges $q^{(a)}{}_{r}$,
$a=1,2,3,4,$ and $5$ (top to bottom) of the amino acid residue types (shown by
their standard one-letter codes) are displayed for KH (red) and Mpipi (green).
For KH, $\varepsilon^{\rm KH}_{r,r^{\prime}}$ values [horizontal variable in
(a)] are taken from Table S3 of ref. (40) (corresponding to the ‘‘KH-D’’ model
in this reference). For Mpipi, $\varepsilon^{\rm Mpipi}_{r,r^{\prime}}$ values
[horizontal variable in (b)] are the contact energies given in Supplementary
Table 11 of ref. (111).

Fig. 1 shows significant deviations from $\varepsilon^{\rm KH}_{r,r^{\prime}}$
or $\varepsilon^{\rm Mpipi}_{r,r^{\prime}}$ when only one eigenvalue is used
($\tilde{a}=1$). At the same time, it indicates that eigenvalue decomposition
is reasonably accurate for $\tilde{a}\geq 2$. Interestingly, for
$\tilde{a}=1$, i.e., when only a single nonelectric charge type is used to
describe the interactions, large deviations are observed only for repulsive
($>0$) interactions in KH but observed for strongly attractive (more negative)
interactions as well as a subset of the less attractive interactions in Mpipi.
When these two pairwise interaction schemes are compared, the nonelectric
charges ($q^{(1)}{}_{r}$s) of the large nonpolar residues---valine,
isoleucine, leucine, and methionine---associated with the energy matrix’s most
dominant ($a=1$) eigenvalue is much higher in KH than in Mpipi, reflecting the
weaker interactions among these nonpolar residues ascribed by Mpipi than by
KH. But the corresponding nonelectric charges for the aromatic residues
phenylalanine, tyrosine, and tryptophan are higher in Mpipi than in KH,
reflecting the relatively more favorable interactions enjoyed by these
residues in Mpipi than in KH (see below). As previously recognized,HaoLi1997 ;
chan99 ; chanMIT at least for KH, the nonelectric charge associated with the
most dominant eigenvalue corresponds to a measure of hydrophobicity, with
polar residues having smaller $q^{(1)}{}_{r}$s and large nonpolar residues
having larger $q^{(1)}{}_{r}$s (red bars in the top panel of Fig.1c). While
all $a=1$ nonelectric charges are positive, $q^{(a)}{}_{r}$s can take positive
and negative values for $a\geq 2$. Intuitively, for KH, the nonelectric charge
for $a=2$ corresponds to a measure of polarity, with most polar residues
having $q^{(2)}{}_{r}>0$ and most nonpolar residues having $q^{(2)}{}_{r}<0$;
and the nonelectric charge for $a=3$ appears to correlate strongly with
electric charge, with arginine and lysine having the largest (most positive)
$q^{(3)}{}_{r}$s whereas aspartic and glutamic acids having the most negative
$q^{(3)}{}_{r}$s (second and third panels from the top of Fig.1c). These
observations are tentalizing, but it should be emphasized that the nonelectric
charges are basically mathematical constructs. As such, they do not always
lend themselves to simple interpretations in terms of physical interaction
types. A case in point is that the corresponding trend of $q^{(2)}{}_{r}$ and
$q^{(3)}{}_{r}$ values in Mpipi is less straightforward to interpret in terms
of a particular physical characteristics of the amino acid residues.

Besides KH and Mpipi, we consider three interaction matrices constructed from
hydrophobicity scales that have recently been used for simulations of single-
IDP conformations and biomolecular LLPS. We refer to these interaction
matrices as (i) $\varepsilon^{\rm HPS}_{r,r^{\prime}}$, which is derived in
ref. (40) from an OPLS (optimized potentials for liquid simulations)
forcefield-based hydrophobicity scale,Rossky (ii) $\varepsilon^{\rm
Urry}_{r,r^{\prime}}$, which is derived in ref. (106) from a prior analysis by
Urry et al. of experimental data on heat-induced conformational compaction of
host-guest polypentapeptides,Urry and (iii) $\varepsilon^{\rm
FB}_{r,r^{\prime}}$, which is derived in ref. (107) by optimizing agreement
between coarse-grained molecular dynamics simulations based on putative
hydrophobicity scales with experimental data on the radii of gyration of
single-chain IDP conformational ensembles. These interaction matrices take the
general form

$\varepsilon^{X}_{r,r^{\prime}}=\frac{\tilde{\lambda}_{r}+\tilde{\lambda}_{r^{\prime}}}{2}+\Delta,\quad
X={\rm HPS,Urry,FB},$ (3)

where $\tilde{\lambda}_{r}$ are twenty residue-specific
hydrophobicity/hydropathy values [not to be confused with the eigenvalues
$\lambda_{a}$ in Eq. (1)] and $\Delta$ is an overall shift. Matrices in the
form of Eq. (3) have at most two non-zero eigenvalues, given by
$\lambda_{1,2}=\frac{1}{2}\Bigl{[}\sum_{r}\tilde{\lambda}_{r}+\Delta\pm\sqrt{20\sum_{r}(\tilde{\lambda}_{r}+\Delta)^{2}}\;\Bigr{]}$
A visual impression of the similarities and differences among the five
interaction matrices considered in this work are provided by their depiction
in Fig. 2 as heat maps.

Fig. 2: Heat-map representations of the interaction matrices considered in the
present work: (a) $\varepsilon^{\rm KH}_{r,r^{\prime}}$, (b) $\varepsilon^{\rm
Mpipi}_{r,r^{\prime}}$, (c) $\varepsilon^{\rm HPS}_{r,r^{\prime}}$, (d)
$\varepsilon^{\rm Urry}_{r,r^{\prime}}$, and (e) $\varepsilon^{\rm
FB}_{r,r^{\prime}}$. The color scale ranges from red (most attractive) to blue
(least attractive or most repulsive). The HPS, Urry and FB matrices follow
from Eq. (3) where the values of $\tilde{\lambda}_{r}$ are taken from Table S1
in ref. (40) for HPS, Table S2 in ref. (106) for Urry, and Table S7 in ref.
(107) for FB. Overall shift $\Delta=0$ for both HPS and FB, whereas
$\Delta=0.08$ for Urry. Corresponding information for $\varepsilon^{\rm
KH}_{r,r^{\prime}}$ and $\varepsilon^{\rm Mpipi}_{r,r^{\prime}}$ is provided
in the caption of Fig. 1.

It follows from the above discussion that formulation of the interactions
entailed by HPS, Urry, and FB by eigenvalue decomposition in terms of
nonelectric charges are straightforward because it requires at most two types
of nonelectric charges. By comparison, the accuracy of corresponding
formulations using a small number of types of nonelectric charges for
interaction schemes encoded by $20\times 20$ energy matrices that are not
derived from a $20$-value hydrophobicity scale (such as KH and Mpipi) has to
be ascertained. In Fig. 3, the accuracy of the truncated eigenvalue
decompositions of $\varepsilon^{\rm KH}_{r,r^{\prime}}$ and $\varepsilon^{\rm
Mpipi}_{r,r^{\prime}}$ is assessed by considering the quantity

$\Delta\varepsilon_{r,r^{\prime}}=\sum_{a=1}^{\tilde{a}}\lambda_{a}q^{(a)}{}_{r}q^{(a)}{}_{r^{\prime}}-\varepsilon_{r,r^{\prime}},$
(4)

which is the difference between the truncated summation and the original
$\varepsilon_{r,r^{\prime}}$. Results for KH and Mpipi at $\tilde{a}=1,2,$ and
$3$ are shown, respectively, in Fig. 3a, b, and c. The heat maps in Fig. 3
indicate that the KH interaction matrix is very well approximated at
$\tilde{a}=3$ with the only exception of the cysteine--cysteine (C--C) entry,
which is noticeably overestimated at all $\tilde{a}\leq 3$. This exception is
in line with the fact that the covalent C--C disulfide bond is physically
different from all other pairwise residue contacts that are noncovalent in
nature. Disulfide bonds are significantly more important in folded proteins
than in IDPs. This may be a reason why the C--C exception in the accuracy of
eigenvalue decomposition is observed in PDB-derived KH but not in the IDP
physics/bioinformatics-derived Mpipi. For Mpipi, the truncated eigenvalue
decomposition at $\tilde{a}=3$ is reasonably accurate though it overestimates
interactions involving the two residues, aspartic and glutamic acids, that
carry negative electric charges.

Under the next subheading, we will use the eigenvalue decomposition in Eq. (1)
to derive a statistical field theory for protein solutions with spatially
short-range residue-dependent interactions governed by a matrix of contact
energies $\varepsilon_{r,r^{\prime}}$. Relative to the field theory we
consider previously for polyampholytes,Pal2021 ; mimb2022 ; wessen2021 ;
wessen2022 the present extended field theory contains $\tilde{a}$ additional
fields, i.e., one additional field $\varphi_{a}({\bf r})$ for every eigenvalue
$\lambda_{a}$ included, wherein the field $\varphi_{a}({\bf r})$ is conjugate
to the density of nonelectric charges $q^{(a)}{}_{r}$. Under subsequent
subheadings, the new theory is first studied using an approximate analytical
approach we find useful for computationally efficient comparisons of
interaction schemes encoded by different $\varepsilon_{r,r^{\prime}}$s. The
theory is then explored using FTS to take into complete account of field
fluctuations so as to provide, in the context of the present model, full
sequence dependence of biomolecular LLPS and detailed structural information
about residue partitioning inside protein-dense condensates.

Fig. 3: Accuracy of the truncated eigenvalue decomposition of $20\times 20$
interaction matrices that are not derived from 20-value hydrophobicity scales.
The level of accuracy for at $\tilde{a}=1,2,$ and $3$ (a--c) is quantified by
$\Delta\varepsilon_{r,r^{\prime}}$ in Eq. (4). Bottom-left triangles using the
blue--red color scale on the left display $\Delta\varepsilon_{r,r^{\prime}}$
for KH, whereas top-right triangles using the green-orange color scale on the
right display the same quantity for Mpipi.

Statistical Field Theory. Consider a system of volume $V$ with $n_{\rm p}$
protein chains, each consisting of $N$ consecutive amino acid residues of
types $r_{\alpha}$ (with $\alpha=1,\dots,N$), and denote the position of
residue $\alpha$ on chain $i$ ($=1,\dots,n_{\rm p}$) as ${\bf R}_{i,\alpha}$.
The system contains also $n_{+}$ unit-positive and $n_{-}$ unit-negative ions
with positions ${\bf r}_{\pm,i}$, and $n_{\rm w}$ neutral solvent particles
with positions ${\bf r}_{{\rm w},i}$. We let $\rho_{\rm p,\pm,w}\equiv n_{\rm
p,\pm,w}/V$ denote the bulk number densities of the (polymers, ions and
solvent) in the system, and $\rho_{\rm b}=Nn_{\rm p}/V=N\rho_{\rm p}$ the bulk
number density of polymer beads. The canonical partition function of this
system is given by

$Z=\frac{V^{n_{\rm p}+n_{+}+n_{-}+n_{\rm w}}}{n_{\rm p}!n_{+}!n_{-}!n_{\rm
w}!}\int\\{\mathrm{d}{\bf R}\\}\\{\mathrm{d}{\bf r}_{+}\\}\\{\mathrm{d}{\bf
r}_{-}\\}\\{\mathrm{d}{\bf r}_{\rm w}\\}\,\mathrm{e}^{-\hat{H}},$ (5)

where

$\\{\mathrm{d}{\bf R}\\}\equiv\prod_{i=1}^{n_{\rm
p}}\prod_{\alpha=1}^{N}\mathrm{d}{\bf R}_{i,\alpha}\,,\quad\\{\mathrm{d}{\bf
r}_{\pm}\\}\equiv\prod_{i=1}^{n_{\pm}}\mathrm{d}{\bf
r}_{\pm,i}\,,\quad\\{\mathrm{d}{\bf r}_{\rm w}\\}\equiv\prod_{i=1}^{n_{\rm
w}}\mathrm{d}{\bf r}_{{\rm w},i}\,$ (6)

are position integration measures for polymer beads, ions and solvents,
respectively. The microscopic Hamiltonian in units of $k_{\rm B}T$, where
$k_{\rm B}$ is Boltzmann’s constant and $T$ is absolute temperature, written
as

$\hat{H}=\hat{H}_{0}+\hat{H}_{\rm c}+\hat{H}_{\rm e}+\hat{H}_{\rm h}\;,$ (7)

contains terms accounting for chain connectivity ($\hat{H}_{0}$), soft
compressibility ($\hat{H}_{\rm c}$), electrostatic interactions ($\hat{H}_{\rm
e}$) and spatially short-range hydrophobic-like interactions---which may
include van der Waals, $\pi$-related, and other forms of short-spatial-range
interactions ($\hat{H}_{\rm h}$). The chain connectivity term

$\hat{H}_{0}=\frac{3}{2b^{2}}\sum_{i=1}^{n_{\rm p}}\sum_{\alpha=1}^{N-1}({\bf
R}_{i,\alpha+1}-{\bf R}_{i,\alpha})^{2}$ (8)

is controlled by the segment length $b$ which we set to the trans Cα\--Cα
virtual bond length, i.e., $b=3.8$ Å, throughout this work. The soft-
compressibility term

$\hat{H}_{\rm c}=\frac{1}{2\gamma}\int\mathrm{d}{\bf r}\left[\hat{\rho}({\bf
r})-\rho_{0}\right]^{2}$ (9)

penalizes deviations of the total density

$\hat{\rho}({\bf r})=v_{\rm b}\hat{\rho}_{\rm b}({\bf
r})+v_{+}\hat{\rho}_{+}({\bf r})+v_{-}\hat{\rho}_{-}({\bf r})+v_{\rm
w}\hat{\rho}_{\rm w}({\bf r})$ (10)

from a reference density $\rho_{0}$, and thus serves to capture an essential
aspect of excluded volume effects. In the above equation, $\hat{\rho}_{i}({\bf
r})$ with $i={\rm b},+,-,{\rm w}$ (referring to polymer bead, cation, anion,
and solvent particle, respectively) is the number density of species $i$ at
position ${\bf r}$ and the $v_{i}$s are factors that we use to model the
relative volume of the different particles. The number densities are

$\hat{\rho}_{\rm b}({\bf r})=\sum_{i=1}^{n_{\rm
p}}\sum_{\alpha=1}^{N}\Gamma({\bf r}-{\bf
R}_{i,\alpha})\quad\mbox{and}\quad\hat{\rho}_{i}({\bf
r})=\sum_{j=1}^{n_{i}}\Gamma({\bf r}-{\bf r}_{i,j}),\quad i=+,-,{\rm w},$ (11)

where individual particles are modeled as Gaussian distributions $\Gamma({\bf
r})=\mathrm{e}^{-{\bf r}^{2}/2a_{\rm s}^{2}}/(2\pi a_{\rm s}^{2})^{3/2}$,
which serve to regularize infinities arising from particle contact
interactions and self interactions.Wang2010 ; Riggleman2012 We refer to
$\gamma$\---which has units of density---in Eq. (9) as the compressibility,
and note that in the limit of $\gamma\rightarrow 0$, we obtain the standard
incompressibility condition $\mathrm{e}^{-\hat{H}_{\rm
c}}\propto\delta[\hat{\rho}({\bf r})-\rho_{0}]$ where the functional
$\delta$-function enforces $\hat{\rho}({\bf r})=\rho_{0}$ at all spatial
positions ${\bf r}$.

Electrostatic interactions are included in $\hat{H}_{\rm e}$ through the
standard pairwise Coulomb potential,

$\hat{H}_{\rm e}=\frac{l_{\rm B}}{2}\int\mathrm{d}{\bf r}\int\mathrm{d}{\bf
r}^{\prime}\,\frac{\hat{c}({\bf r})\hat{c}({\bf r}^{\prime})}{|{\bf r}-{\bf
r}^{\prime}|}\,,$ (12)

where

$\hat{c}({\bf r})=\sum_{i=1}^{n_{\rm
p}}\sum_{\alpha=1}^{N}\sigma_{i,\alpha}\Gamma({\bf r}-{\bf
R}_{i,\alpha})+\hat{\rho}_{+}({\bf r})-\hat{\rho}_{-}({\bf r}),$ (13)

is the density of electric charge at ${\bf r}$ and $\sigma_{i,\alpha}$ is the
electric charge of the $\alpha$th residue on the $i$th chain in units of the
proton charge $e$. Because the present study is limited to systems of chains
with identical sequence, i.e., $\sigma_{i,\alpha}=\sigma_{\alpha}$ is
independent of $i$, the subscript $i$ will be dropped from $\sigma$ hereafter.
The strength of electrostatic interactions is controlled by the Bjerrum length
$l_{\rm B}=e^{2}/4\pi\epsilon_{0}\epsilon_{\rm r}k_{\rm B}T$ that involves
vacuum permittivity $\epsilon_{0}$ and a constant background relative
permittivity $\epsilon_{\rm r}$. In this work, we neglect any concentration
dependence of $\epsilon_{\rm r}$ as this has been shown to amount to only
small to moderate effects on LLPS propensity.wessen2021 Nonetheless, if
desired, it would be straightforward to implement a concentration dependent
permittivity in our model through the approach described in ref. (89).

The final term in the Hamiltonian in Eq. (7), $\hat{H}_{\rm h}$, contains the
pairwise residue-specific spatially short-range interactions, which we express
as

$\hat{H}_{\rm h}=\frac{1}{2}\sum_{i,j=1}^{n_{\rm
p}}\sum_{\alpha,\beta=1}^{N}\varepsilon_{r_{\alpha},r_{\beta}}\,\Gamma^{2}\star
V_{\rm h}(|{\bf R}_{i,\alpha}-{\bf R}_{j,\beta}|)$ (14)

where

$\Gamma^{2}\star V_{\rm h}(|{\bf R}_{i,\alpha}-{\bf
R}_{j,\beta}|)\equiv\int\mathrm{d}{\bf r}\int\mathrm{d}{\bf
r}^{\prime}\,\Gamma({\bf r}-{\bf R}_{i,\alpha})\Gamma({\bf r}^{\prime}-{\bf
R}_{j,\beta})V_{\rm h}(|{\bf r}-{\bf r}^{\prime}|)$ (15)

accounts for the residue-residue interaction in the presence of Gaussian
smearing (note that $\Gamma^{2}\star V_{\rm h}\rightarrow V_{\rm h}$ as
$a_{\rm s}\rightarrow 0$). Leaving the functional form of the interaction
potential $V_{\rm h}(|{\bf r}|)$ (in units of $k_{\rm B}T$) unspecified for
the moment, we have in the above expression assumed that all residue
dependence are accounted for by the overall multiplicative factors
$\varepsilon_{r_{\alpha},r_{\beta}}$ from a given interaction matrix. The
spectral form of $\varepsilon_{r,r^{\prime}}$ in Eq. (1) can be used to write
$\hat{H}_{\rm h}$ in Eq. (14) as

$\hat{H}_{\rm h}=\sum_{a=1}^{20}\frac{1}{2}\int\mathrm{d}{\bf
r}\int\mathrm{d}{\bf r}^{\prime}\,V_{\rm h}(|{\bf r}-{\bf
r}^{\prime}|)\lambda_{a}\hat{h}_{a}({\bf r})\hat{h}_{a}({\bf r}^{\prime}),$
(16)

where

$\hat{h}_{a}({\bf r})=\sum_{i=1}^{n_{\rm
p}}\sum_{\alpha=1}^{N}q^{(a)}{}_{r_{\alpha}}\Gamma({\bf r}-{\bf
R}_{i,\alpha})$ (17)

is the density of nonelectric charge for the $a$th eigenvalue $\lambda_{a}$ at
position ${\bf r}$.

The Hamiltonian terms $\hat{H}_{\rm c,e,h}$ for non-bonded interactions are
quadratic in densities and can thus be decoupled using standard Hubbard-
Stratonovich transformations,Fredrickson2006 turning the partition function in
Eq. (5) into that of a statistical field theory. Special care needs to be
given to $\hat{H}_{\rm h}$ wherein terms with $\lambda_{a}>0$ and
$\lambda_{a}<0$ need to be linearized with fields integrated on the imaginary
and real axis, respectively. To keep track of the sign of eigenvalue
$\lambda_{a}$, we introduce the variable $\xi_{a}$, defined as

$\xi_{a}=\left\\{\begin{matrix}-\mathrm{i}&,&\lambda_{a}<0\\\
1&,&\lambda_{a}>0\end{matrix}\right.,$ (18)

where $\mathrm{i}^{2}=-1$. The partition function $Z$ in Eq. (5) can then be
shown to be equivalent to a statistical field theory with partition function

$Z=\frac{V^{n_{\rm p}+n_{+}+n_{-}+n_{\rm w}}}{n_{\rm p}!n_{+}!n_{-}!n_{\rm
w}!}\int{\cal D}\eta\int{\cal D}\psi\left(\prod_{a=1}^{20}\int{\cal
D}\varphi_{a}\right)\mathrm{e}^{-H[\eta,\psi,\\{\varphi_{a}\\}]},$ (19)

with field Hamiltonian (in units of $k_{\rm B}T$)

$H=-\sum_{i={\rm p,+,-,w}}n_{i}\ln
Q_{i}[\breve{\eta},\breve{\psi},\\{\breve{\varphi}_{a}\\}]+\int\mathrm{d}{\bf
r}\left[-\mathrm{i}\rho_{0}\eta+\frac{\gamma\eta^{2}}{2}+\frac{(\bm{\nabla}\psi)^{2}}{8\pi
l_{\rm B}}+\sum_{a=1}^{20}\frac{\varphi_{a}V_{\rm
h}{}^{-1}\varphi_{a}}{2|\lambda_{a}|}\right].$ (20)

In the field theory representation of the system, $\eta({\bf r})$ and
$\psi({\bf r})$ are conjugate fields to the density $\hat{\rho}({\bf r})$ and
the electric charge density $\hat{c}({\bf r})$, respectively, while
$\varphi_{a}({\bf r})$ is conjugate to the nonelectric charge density
$\hat{h}_{a}({\bf r})$ associated with eigenvalue $\lambda_{a}$ in the
spectral decomposition of the interaction matrix. The single-molecule
partition functions $Q_{\rm
p,\pm,w}[\breve{\eta},\breve{\psi},\\{\breve{\varphi}_{a}\\}]$ depend on
smeared fields $\breve{\phi}({\bf r})=\Gamma\star\phi({\bf
r})\equiv\int\mathrm{d}{\bf r}^{\prime}\Gamma({\bf r}-{\bf
r}^{\prime})\phi({\bf r}^{\prime})$ (for $\phi=\eta,\psi,\varphi_{a}$) as

$\displaystyle Q_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{V}\int\mathrm{d}{\bf
r}\exp\left[-\mathrm{i}(v_{\pm}\breve{\eta}({\bf r})\pm\breve{\psi}({\bf
r}))\right]\,,$ (21) $\displaystyle Q_{\rm w}$ $\displaystyle=$
$\displaystyle\frac{1}{V}\int\mathrm{d}{\bf r}\exp\left[-\mathrm{i}v_{\rm
w}\breve{\eta}({\bf r})\right]\,$ (22)

for ions and solvents, and

$Q_{\rm p}=\frac{1}{\mathcal{N}}\left(\prod_{\alpha=1}^{N}\int\mathrm{d}{\bf
R}_{\alpha}\right)\exp\left[-\frac{3}{2b^{2}}\sum_{\alpha=1}^{N-1}({\bf
R}_{\alpha+1}-{\bf R}_{\alpha})^{2}-\sum_{\alpha=1}^{N}W_{\alpha}({\bf
R}_{\alpha})\right]\,$ (23)

for polymer chains, where $\mathcal{N}\equiv V(2\pi b^{2}/3)^{3(N-1)/2}$
normalizes $Q_{\rm p}$ to unity at zero field values, and

$W_{\alpha}({\bf r})=\mathrm{i}v_{\rm b}\breve{\eta}({\bf
r})+\mathrm{i}\sigma_{\alpha}\breve{\psi}({\bf
r})+\sum_{a=1}^{20}\mathrm{i}\,\xi_{a}q^{(a)}{}_{r_{\alpha}}\breve{\varphi}_{a}({\bf
r})\,.$ (24)

In Eq. (20), $V_{\rm h}{}^{-1}$ represents a differential operator satisfying

$V_{\rm h}{}^{-1}V_{\rm h}(|{\bf r}|)=\delta({\bf r})$ (25)

and thus $V_{\rm h}{}^{-1}$ depends on the functional form of $V_{\rm h}(|{\bf
r}|)$.

Analytical Calculation of the Leading Order Free Energy. We next proceed to an
approximate evaluation of the statistical field theory partition function $Z$
in Eq. (19). In the analytical approach taken here, we only account for the
leading-order effects from each field around the homogeneous saddle solution.
For our model, this amounts to the mean field theory (MFT) solution for the
fields $\eta$ and $\varphi_{a}$, while electrostatic effects first show up at
the one-loop level in $\psi$ which are accounted for by RPA.

The spatially homogeneous MFT solution is obtained by first letting the fields
$\eta$, $\psi$ and $\varphi_{a}$ in the expression for the field Hamiltonian
$H$ in Eq. (20) be position-independent (hence smeared fields
$\breve{\phi}=\phi$ for $\phi=\eta,\psi,\varphi_{a}$ are also position-
independent) and then requiring that the first derivatives of the resulting
$H$ with respect to the fields are zero. With position-independent fields, the
single molecule partition functions $Q_{i}$, defined in Eqs. (21)---(23),
become

$\displaystyle Q_{\rm w}$ $\displaystyle=\mathrm{e}^{-\mathrm{i}v_{\rm
w}\eta},$ (26) $\displaystyle Q_{\pm}$
$\displaystyle=\mathrm{e}^{-\mathrm{i}v_{\pm}\eta\mp\mathrm{i}\psi},$
$\displaystyle Q_{\rm p}$ $\displaystyle=\exp\left[-\mathrm{i}Nv_{\rm
b}\eta-\mathrm{i}\sigma_{\rm
p}\psi-\mathrm{i}\sum_{a=1}^{20}\xi_{a}q^{(a)}\varphi_{a}\right],$

where $\sigma_{\rm p}=\sum_{\alpha=1}^{N}\sigma_{\alpha}$ is the net charge of
a polymer and $q^{(a)}=\sum_{\alpha=1}^{N}q^{(a)}{}_{r_{\alpha}}$ is the total
nonelectric charge associated with the $a$th eigenvalue of one polymer chain.
The above relations for $Q_{i}$ lead to the following expression for $H$ with
homogeneous fields,

$\displaystyle H=$ $\displaystyle\mathrm{i}(n_{\rm p}Nn_{\rm
b}+n_{+}v_{+}+n_{-}v_{-}+n_{\rm w}v_{\rm w})\eta+\mathrm{i}(n_{\rm
p}\sigma_{\rm p}+n_{+}-n_{-})\psi+\mathrm{i}n_{\rm
p}\sum_{a=1}^{20}\xi_{a}q^{(a)}\varphi_{a}$ (27)
$\displaystyle+V\left[-\mathrm{i}\rho_{0}\eta+\frac{\gamma\eta^{2}}{2}+\sum_{a=1}^{20}\frac{\varphi_{a}^{2}}{2|\lambda_{a}|\hat{V}_{\rm
h}(\bm{0})}\right].$

In Eq. (27), $\hat{V}_{\rm h}(\bm{0})$ is the ${\bf k}=\bm{0}$ component of
the Fourier transform $\hat{V}_{\rm h}({\bf k})=\int\mathrm{d}{\bf
r}\,\mathrm{e}^{-\mathrm{i}{\bf r}\cdot{\bf k}}V_{\rm h}(|{\bf r}|)$ of the
spatially short-range potential $V_{\rm h}(|{\bf r}|)$. This factor can be
derived by considering the Fourier representation of $\delta({\bf
r})=(2\pi)^{-3}\int\mathrm{d}{\bf k}\,\exp(\mathrm{i}{\bf k}\cdot{\bf r})$ in
the defining Eq. (25) for $V_{\rm h}^{-1}$, which leads to $V_{\rm
h}^{-1}\exp(\mathrm{i}{\bf k}\cdot{\bf r})=\exp(\mathrm{i}{\bf k}\cdot{\bf
r})/\hat{V}_{\rm h}({\bf k})$. Then, by noting that the Fourier transform of
the position-independent $\varphi_{a}$ is proportional to $\delta({\bf k})$,
it is straightforward to show that the $\int\mathrm{d}{\bf r}$ integration
with position-independent $\varphi_{a}$ produces the $\hat{V}_{\rm h}(\bm{0})$
term.

The coefficient of $\mathrm{i}\psi$ in Eq. (27) is the total electric charge
in the system which must be set to zero, $n_{\rm p}\sigma_{\rm
p}+n_{+}-n_{-}=0$. The charge-conjugate field $\psi$ thus vanishes in Eq. (27)
and does not contribute to the MFT solution. For the other fields, solving
$\partial H/\partial\eta=\partial H/\partial\varphi_{a}=0$ for $\bar{\eta}$
and $\bar{\varphi}_{a}$ leads to

$\bar{\eta}=\frac{-\mathrm{i}\Delta\rho}{\gamma}\quad\mbox{and}\quad\bar{\varphi}_{a}=-\mathrm{i}\hat{V}_{\rm
h}(\bm{0})\xi_{a}|\lambda_{a}|q^{(a)}\rho_{\rm p},$ (28)

where $\Delta\rho\equiv\nu_{\rm b}N\rho_{\rm
p}+\nu_{+}\rho_{+}+\nu_{-}\rho_{-}+\nu_{\rm w}\rho_{\rm w}-\rho_{0}$ is the
deviation of the total bulk density from $\rho_{0}$, and
$q^{(a)}\equiv\sum_{\alpha=1}^{N}q^{(a)}{}_{r_{\alpha}}$ is the total
nonelectric charge associated with the $a$th eigenvalue $\lambda_{a}$ of one
polymer chain. The same expressions for $\bar{\eta}$ and $\bar{\varphi}_{a}$
[Eq. (28)] may alternatively be obtained by requiring vanishing first
functional derivatives of $H$ with respect to $\eta({\bf r})$ and
$\varphi_{a}({\bf r})$ and then seeking position-independent solutions
$\eta({\bf r})\rightarrow\bar{\eta}$, $\psi({\bf r})\rightarrow\bar{\psi}$ and
$\varphi_{a}({\bf r})\rightarrow\bar{\varphi}_{a}$ of the resulting $\delta
H/\delta\eta({\bf r})=\delta H/\delta\psi({\bf r})=\delta
H/\delta\varphi_{a}({\bf r})=0$ conditions, as will be discussed below under
the next subheading.

Plugging the MFT solution Eq. (28) back into Eq. (27) leads to the field
Hamiltonian per unit volume $\bar{h}=H/V$, given by

$\bar{h}=\frac{{\Delta\rho}^{2}}{2\gamma}+\frac{\hat{V}_{\rm
h}(\bm{0})}{2}\left[\sum_{a=1}^{20}\lambda_{a}({q^{(a)}})^{2}\right]\rho_{\rm
p}^{2}\;.$ (29)

In the approximate analytical approach taken in this work, we set $\eta({\bf
r})$ and $\varphi_{a}({\bf r})$ to their MFT values $\bar{\eta}$ and
$\bar{\varphi}_{a}$, while a spatially varying $\psi({\bf r})$ is kept to
$\mathcal{O}(\psi^{2})$ in $H$. The field Hamiltonian in this approximation is
therefore $H\approx V\bar{h}+H_{\rm RPA}[\psi]$ where the
$\mathcal{O}(\psi^{2})$ terms are contained in

$H_{\rm RPA}=\int\frac{\mathrm{d}{\bf
k}}{(2\pi)^{3}}\frac{1}{2}\hat{\psi}(-{\bf
k})\left[\hat{\Gamma}^{2}\left(N\rho_{\rm p}g_{\rm cc}(|{\bf
k}|)+\rho_{+}+\rho_{-}\right)+\frac{{\bf k}^{2}}{4\pi l_{\rm
B}}\right]\hat{\psi}({\bf k}).$ (30)

The factor $g_{\rm cc}(k)$, given bywessen2021

$g_{\rm
cc}(k)=\frac{1}{N}\sum_{\alpha,\beta=1}^{N}\sigma_{\alpha}\sigma_{\beta}\mathrm{e}^{-|\alpha-\beta|b^{2}k^{2}/6},$
(31)

comes from the quadratic expansion of the single polymer partition function
$Q_{\rm p}$. In our approximation, the functional integrals in the partition
function $Z$ amounts to

$\int{\cal D}\eta\int{\cal D}\psi\left(\prod_{a=1}^{20}\int{\cal
D}\varphi_{a}\right)\mathrm{e}^{-H[\eta,\psi,\\{\varphi_{a}\\}]}\approx\mathrm{e}^{-V\bar{h}}\int{\cal
D}\psi\,\mathrm{e}^{-H_{\rm RPA}}$ (32)

where $V\bar{h}$ is the field Hamiltonian evaluated at homogeneous field
values in accordance with Eq. (29).

To simplify the formulation, we now take the $\gamma\rightarrow 0$ limit such
that solvent density is no longer an independent component of the system but
is instead determined by

$\rho_{\rm w}=\frac{1}{v_{\rm w}}\left(\rho_{0}-v_{\rm b}N\rho_{\rm
p}-v_{+}\rho_{+}-v_{-}\rho_{-}\right).$ (33)

In other words, $\Delta\rho=0$ and, equivalently, $\rho_{0}=v_{\rm
b}N\rho_{\rm p}+v_{+}\rho_{+}+v_{-}\rho_{-}+v_{\rm w}\rho_{\rm w}$. The
remaining system components $\\{\rho_{\rm p},\rho_{+},\rho_{-}\\}$ are further
constrained by electric charge neutrality of the system as a whole, i.e.
$\sigma_{\rm p}\rho_{\rm p}+\rho_{+}-\rho_{-}=0$, where $\sigma_{\rm p}$ is
the net electric charge of a polymer chain. To make this condition manifest,
we introduce the overall salt density $\rho_{\rm s}$, which is the overall
number density of cation-anion pairs. The ion type in excess is referred to as
counterions. For net-neutral or net-positive chains, $\sigma_{\rm p}\geq 0$,
we have $\rho_{+}=\rho_{\rm s}$ and $\rho_{-}=\sigma_{\rm p}\rho_{\rm
p}+\rho_{\rm s}$; for net-negative chains, $\sigma_{\rm p}<0$, we have we have
$\rho_{-}=\rho_{\rm s}$ and $\rho_{+}=-\sigma_{\rm p}\rho_{\rm p}+\rho_{\rm
s}$. Taking all the above considerations together and with $k\equiv|{\bf k}|$,
the MFT/RPA free energy per volume in units of $k_{\rm B}T$, $f(\rho_{\rm
p},\rho_{\rm s})\equiv-(\ln Z)/V$, becomes

$f(\rho_{\rm p},\rho_{\rm s})=-s(\rho_{\rm p},\rho_{\rm s})-\rho_{0}\chi_{\rm
eff}\phi_{\rm
p}^{2}+\frac{1}{4\pi}\int_{0}^{\infty}\mathrm{d}k\,k^{2}\ln\left[1+\frac{4\pi
l_{\rm B}}{k^{2}}\hat{\Gamma}^{2}\left(2\rho_{\rm s}+|\sigma_{\rm p}|\rho_{\rm
p}+g_{\rm cc}(k)\rho_{\rm p}\right)\right]$ (34)

after performing the $\int{\cal D}\psi$ integrationwessen2021 in Eq. (32),
wherein the term

$-s(\rho_{\rm p},\rho_{\rm s})=\sum_{i={\rm p,+,-,w}}\rho_{i}\ln\rho_{i}\,$
(35)

accounts for translation entropy and follows from applying Stirling’s
approximation $\ln n!\approx n\ln n-n$ to the factorial prefactors in Eq. (19)
and neglect terms linear in $\rho_{i}$ arising from the nonlogarithmic $-n$
part of the Stirling approximation because terms linear in $\rho_{i}$ in the
free energy have no effect on phase separation properties.linJML In the second
term on the right hand side of Eq. (34), $\phi_{\rm p}\equiv v_{\rm
b}N\rho_{\rm p}/\rho_{0}$ is the volume fraction occupied by polymer beads
(and thus the volume fraction of the polymers themselves) and $\chi_{\rm eff}$
is a dimensionless effective Flory-Huggins $\chi$-parameter that originates
from the spatially short-range interactions [cf. second term for $\bar{h}$ in
Eq. (29)]:

$\chi_{\rm eff}=-\frac{\rho_{0}\hat{V}_{\rm h}(\bm{0})}{2v_{\rm
b}^{2}N^{2}}\sum_{a=1}^{20}\lambda_{a}({q^{(a)}})^{2}\,.$ (36)

Using Eq. (1), the summation over $a$ in the above expression for $\chi_{\rm
eff}$ can be recast,

$\sum_{a=1}^{20}\lambda_{a}({q^{(a)}})^{2}=\sum_{\alpha=1}^{N}\sum_{\beta=1}^{N}\varepsilon_{r_{\alpha},r_{\beta}},$
(37)

to show that it is a summation over all possible pairwise residue-residue
interaction energies for a pair of heteropolymer chains with the given
sequence of residues. As such, the present mean-field treatment of spatially
short-range interactions is akin to the random-mixing or Bragg-Williamsbragg-
williams approximation used in Flory-Huggins theories of proteinlike
heteropolymers.dill1985 As mentioned above, the third term in Eq. (34)
corresponds to the standard RPA term accounting for Gaussian fluctuations in
electric charge density. This term is obtained from functional integral over
$\psi$ in Eq. (32) followed by subtraction of the free-energy contribution at
$\rho_{\rm p}=\rho_{\rm s}=0$, a subtraction that has no effect on phase
separation properties. Given the free energy $f(\rho_{\rm p},\rho_{\rm s})$ in
Eq. (34), phase diagrams may be constructed, for example, by matching the
polymer and salt chemical potentials $\mu_{\rm p,s}=\partial
f/\partial\rho_{\rm p,s}$ and osmotic pressure $\Pi=\mu_{\rm p}\rho_{\rm
p}+\mu_{\rm s}\rho_{\rm s}-f$ (ref. (78)).

Field-Theoretic Simulations (FTS). One shortcoming of the above analytical
approximation is that all effects of the spatially short-range interactions
are condensed into the value of $\chi_{\rm eff}$ which depends only on overall
residue composition of the polymer but not the specific sequential arrangement
of the residues. To tackle such sequence-specific effects, one needs to
account for higher-order fluctuations in the partition function in Eq. (19).
To this end, we next consider ways to study systems described by Eq. (19)
using field FTS because of its ability to afford, in principle, a full account
for field fluctuations.Fredrickson2002

In FTS, each field (denoted generically as $\phi$) is analytically continued
into its complex plane, and approximated by a set of discrete variables
defined on the sites of a cubic $M^{3}$ lattice with periodic boundary
conditions and lattice spacing $\Delta x$. The fields evolve in a fictitious
complex-Langevin (CL) time $t$ according to

$\frac{\partial\phi({\bf r},t)}{\partial t}=-\frac{\delta H}{\delta\phi({\bf
r},t)}+\tilde{\eta}_{\phi}({\bf r},t)\quad,\quad\phi=\eta,\psi,\varphi_{a},$
(38)

where $\tilde{\eta}_{\phi}({\bf r},t)$ represents real-valued Gaussian noise
with zero mean, i.e., $\langle\tilde{\eta}_{\phi}({\bf
r},t)\tilde{\eta}_{\phi^{\prime}}({\bf
r}^{\prime},t^{\prime})\rangle=2\delta_{\phi,\phi^{\prime}}\delta(t-t^{\prime})\delta({\bf
r}-{\bf r}^{\prime})$. Thermal averages of the real-space system are then be
computed as asymptotic CL time averages of field operators in the fictitious-
time system that are constructed to correspond to the thermodynamic observable
of interest.Parisi1983 ; Klauder1983 . This approach has its origin in the
development, beginning in the 1980s, of stochastic quantizationParisiWu1981 as
a method for studying quantum field theories and their
regularizationHSCMartyGhost ; HSCMartyGravity ; rumpf and has since been
applied extensively to study properties of polymer solutions.Fredrickson2006

In the present model, the functional derivatives of the field Hamiltonian $H$
in Eq. (20) are

$\displaystyle\frac{\delta H}{\delta\eta({\bf r})}$
$\displaystyle=\mathrm{i}\sum_{i={\rm,b,\pm,w}}v_{i}\tilde{\rho}_{i}({\bf
r})-\mathrm{i}\rho_{0}+\gamma\eta({\bf r}),$ (39) $\displaystyle\frac{\delta
H}{\delta\psi({\bf r})}$ $\displaystyle=\mathrm{i}\tilde{c}({\bf
r})-\frac{\bm{\nabla}^{2}\psi({\bf r})}{4\pi l_{\rm B}},$
$\displaystyle\frac{\delta H}{\delta\varphi_{a}({\bf r})}$
$\displaystyle=\mathrm{i}\xi_{a}\tilde{h}_{a}({\bf
r})+\frac{1}{|\lambda_{a}|}V_{\rm h}^{-1}\varphi_{a}({\bf r})\;,$

where $\tilde{\rho}_{i}({\bf r})$ in the first relation here in Eq. (39) is a
field operator corresponding to the density of component $i$, and is given,
respectively, by

$\tilde{\rho}_{\pm}({\bf
r})=\Gamma\star\frac{\rho_{\pm}}{Q_{\pm}}\mathrm{e}^{-\mathrm{i}v_{\pm}\breve{\eta}({\bf
r})\mp\mathrm{i}\breve{\psi}({\bf r})}\quad\mbox{and}\quad\tilde{\rho}_{\rm
w}({\bf r})=\Gamma\star\frac{\rho_{\rm w}}{Q_{\rm
w}}\mathrm{e}^{-\mathrm{i}v_{\rm w}\breve{\eta}({\bf r})}$ (40)

for ions and solvents, whereas the polymer bead density operator,
$\tilde{\rho}_{\rm b}({\bf r})$, is computed through forward and backward
chain propagators $q_{\rm F,B}({\bf r},\alpha)$ which are constructed by
applying the following relations iteratively:Pal2021 ; mimb2022 ; wessen2021

$\displaystyle q_{\rm F}({\bf r},\alpha+1)=\mathrm{e}^{-W_{\alpha+1}({\bf
r})}\,\Phi\star q_{\rm F}({\bf r},\alpha),$ (41) $\displaystyle q_{\rm B}({\bf
r},\alpha-1)=\mathrm{e}^{-W_{\alpha-1}({\bf r})}\,\Phi\star q_{\rm B}({\bf
r},\alpha),$

where $\Phi({\bf r})=\mathrm{e}^{-3{\bf r}^{2}/2b^{2}}/(2\pi b^{2}/3)^{3/2}$,
and initiating from $q_{\rm F}({\bf r},1)=\mathrm{e}^{-W_{1}({\bf r})}$ and
$q_{\rm B}({\bf r},N)=\mathrm{e}^{-W_{N}({\bf r})}$. Given $q_{\rm F,B}$, the
polymer bead density operator can be computed as

$\tilde{\rho}_{\rm b}({\bf r})=\Gamma\star\frac{\rho_{\rm p}}{Q_{\rm
p}}\sum_{\alpha=1}^{N}\mathrm{e}^{W_{\alpha}({\bf r})}q_{\rm F}({\bf
r},\alpha)q_{\rm B}({\bf r},\alpha),$ (42)

with

$Q_{\rm p}=\frac{1}{V}\int\mathrm{d}{\bf r}\,q_{\rm F}({\bf r},N).$ (43)

In the second relation in Eq. (39), $\tilde{c}({\bf r})=\tilde{c}_{\rm b}({\bf
r})+\tilde{\rho}_{+}({\bf r})-\tilde{\rho}_{-}({\bf r})$ is a field operator
corresponding to the electric charge density, wherein the polymer bead
contribution $\tilde{c}_{\rm b}({\bf r})$ is given by

$\tilde{c}_{\rm b}({\bf r})=\Gamma\star\frac{\rho_{\rm p}}{Q_{\rm
p}}\sum_{\alpha=1}^{N}\mathrm{e}^{W_{\alpha}({\bf r})}q_{\rm F}({\bf
r},\alpha)q_{\rm B}({\bf r},\alpha)\sigma_{\alpha}.$ (44)

In the third relation in Eq. (39), the field operators for the nonelectric
charge densities associated with the spatially short-range interactions,
$\tilde{h}_{a}({\bf r})$, are similarly given by

$\tilde{h}_{a}({\bf r})=\Gamma\star\frac{\rho_{\rm p}}{Q_{\rm
p}}\sum_{\alpha=1}^{N}\mathrm{e}^{W_{\alpha}({\bf r})}q_{\rm F}({\bf
r},\alpha)q_{\rm B}({\bf r},\alpha)q^{(a)}{}_{r_{\alpha}}\;.$ (45)

It can be readily verified that the position-independent solution to $\delta
H/\delta\phi({\bf r})=0$ for $\phi=\eta,\psi,\varphi_{a}$ (where the
expressions for the functional derivatives are given in Eq. (39)) yields
exactly the MFT solution in Eq. (28) for position-independent $\eta$ and
$\varphi_{a}$ (together with any value of position-independent $\psi$): For
any values of the position-independent fields, the field operators for number
and charge densities become equal to their bulk counter-parts, i.e.
$\tilde{\rho}_{\rm b,\pm,w}({\bf r})=\rho_{\rm b,\pm,w}$, $\tilde{c}({\bf
r})=\rho_{\rm p}\sigma_{\rm p}+\rho_{+}-\rho_{-}=0$ and $\tilde{h}_{a}({\bf
r})=\rho_{\rm p}q^{(a)}$. Substituting this into Eq. (39) and setting the
functional derivatives to zero gives exactly the expressions in Eq. (28).

In contrast to the approximate analytical approach introduced under the
previous subheading, FTS depends on the full analytical form of the spatially
short-range interaction potential $V_{\rm h}(r)$, where $r=|{\bf r}|$, not
only a factor of its spatial integral $\int\mathrm{d}{\bf r}V_{\rm
h}(r)=\hat{V}_{\rm h}(\bm{0})$ that appears as part of an effective Flory-
Huggins $\chi$-parameter, $\chi_{\rm eff}$, in Eq. (36). Several functional
forms for spatially short-range interactions have been used extensively in
recent coarse-grained simulations of IDP LLPS.dignon18 ; suman2 ; Mpipi These
include the common Lennard-Jones (LJ) potential

$V_{\rm LJ}({\cal E},\sigma|r)=4{\cal
E}\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]\;$
(46)

where ${\cal E}$ is the depth of the potential well and $\sigma$ is a length
scale, LJ variations such as the shifted Weeks–Chandler–Andersen (WCA)
potentialWCA for finite-range purely repulsive interactions, and the recently
proposed Wang-Frenkel potentialWangFrenkel with a finite spatial range $r_{\rm
c}$,

$V_{\rm WF}({\cal E},\mu,\nu,r_{\rm c},\sigma|r)={\cal E}\alpha_{\rm
WF}(\mu,\nu,r_{\rm
c},\sigma)\left[\left(\frac{\sigma}{r}\right)^{2\mu}-1\right]\left[\left(\frac{r_{\rm
c}}{r}\right)^{2\mu}-1\right]^{2\nu}\;$ (47)

for $r\leq r_{\rm c}$ and $V_{\rm WF}(\mu,\nu,r_{\rm c},\sigma|r)=0$ for
$r>r_{\rm c}$ and where

$\alpha_{\rm WF}(\mu,\nu,r_{\rm c},\sigma)=2\nu\left(\frac{r_{\rm
c}}{\sigma}\right)^{2\mu}\left\\{\frac{1+2\nu}{2\nu[(r_{\rm
c}/\sigma)^{2\mu}-1]}\right\\}^{2\nu+1}\;,$ (48)

that enjoys several apparent numerical advantages when applied to many-body
simulations, including its smooth decay to zero as $r$ is increased toward the
cutoff distance $r_{\rm c}$. Note that the original notation in ref. (141) is
largely followed in Eqs. (47) and (48). The symbols $\alpha_{\rm WF}$ and
$\mu$ here should not be confused with the polymer bead label
$\alpha=1,2,\dots,N$ and the symbol for chemical potential. The length scale
$\sigma$ in Eqs. (46)--(48) also should not be confused with the symbol for
electric charge in our formulation.

Mathematically, however, these potential functions cannot practically---if at
all possibly---be incorporated into our field-theoretic formalism because the
field Hamiltonian $H$ in Eq. (20) requires the inverse operator, ${V_{\rm
h}}^{-1}$, of the potential function but no inverse operator expressed in
terms of elementary functions is known for $V_{\rm LJ}$ and $V_{\rm WF}$. (For
a promising approach for including general pairwise interaction potentials in
FTS, see ref. (142)). This consideration leads us to the Yukawa
potential,yukawa which has the same mathematical form as the screened Coulomb
potential, because it possesses an inverse expressible in terms of Fourier
transform of elementary functions. Therefore, as a first step in the
exploration of our general theoretical framework, we specialize here to a
Yukawa form for the spatially short-range interactions, viz.,

$V_{\rm h}(r)=\frac{l_{\rm h}\mathrm{e}^{-\kappa r}}{r}\,,$ (49)

where interaction strength is represented by $l_{\rm h}$ and the spatial range
of the interaction is characterized by the reciprocal $\kappa^{-1}$ of the
screening coefficient $\kappa$ (Fig. 4). The Fourier transform of $V_{\rm
h}(r)$ is given by $\hat{V}_{\rm h}({\bf k})=\hat{V}_{\rm h}(k)=4\pi l_{\rm
h}/(k^{2}+\kappa^{2})$, and the inverse of $V_{\rm h}(r)$ as defined by Eq.
(25) is ${V_{\rm h}}^{-1}=(-\bm{\nabla}^{2}+\kappa^{2})/4\pi l_{\rm h}$. In
view of the MFT result in Eq. (36), the overall spatial range-dependent
interaction strength of the Yukawa potential may also be quantified by
$\hat{V}_{\rm h}(\bm{0})=4\pi l_{\rm h}/\kappa^{2}$.

Fig. 4: The Yukawa potential as a model of spatially short-range favorable
interactions. Two examples of $V_{\rm h}(r)$ in Eq. (49) are depicted by the
thick solid black curve ($l_{\rm h}=6\sigma/3.8$, $\kappa=3.8/8\sigma$) and
the thick solid blue curve ($l_{\rm h}=3\sigma/3.8$, $\kappa=3.8/5\sigma$)
with $\sigma=1$; i.e., $\sigma$ is the unit for $r$. The dashed curves (same
color code) provide the corresponding ‘‘smeared’’ version of these potentials
defined as
$(1/2\pi^{2})\int_{0}^{\infty}\mathrm{d}kk^{2}\hat{\Gamma}(k)^{2}\hat{V}(k)\sin(kr)/kr$
where $\hat{\Gamma}(k)=\exp(-a_{\rm s}^{2}k^{2}/2)$ is the Fourier transform
of the function $\Gamma({\bf r})$ in Eq. (11) with $a_{\rm s}=\sigma/\sqrt{6}$
shown as an example here. The interaction between two point particles via this
smeared Yukawa potential is mathematically equivalent to that between two
particles modeled as Gaussian distributions $\Gamma({\bf r})$ via the original
$V_{\rm h}(r)$. Included for comparison are a LJ potential (red curve, Eq.
(46) with ${\cal E}=1$) and a Wang-Frenkel potential (thin blue curve, Eq.
(47) with ${\cal E}=1$, $\mu=\nu=1$, and $r_{\rm c}=3\sigma$). The example
Yukawa potentials are selected to match approximately with either the
potential wells (thick blue curves) or the large-$r$ trends (thick black
curves) of the LJ and Wang-Frenkel potentials.

Structural information of the molecular species in our model system can be
inferred from pair correlation functions,Pal2021 ; wessen2022 defined here as

$G_{p,q}(|{\bf r}-{\bf r}^{\prime}|)=\langle\hat{\rho}_{p}({\bf
r})\hat{\rho}_{q}({\bf
r}^{\prime})\rangle-\delta_{pq}\rho_{q}\frac{\mathrm{e}^{-|{\bf r}-{\bf
r}^{\prime}|^{2}/4a_{\rm s}^{2}}}{(4\pi a_{\rm s}^{2})^{3/2}}\;,\quad\quad
p,q={\rm b,\pm,w}\;,$ (50)

where the second term on the right hand side of the equation for $G_{p,q}$
subtracts the trivial correlation of a particle with itself. Correlation
functions were used to quantify the degree of condensed-state
subcompartmentalization of phase-separated polyampholytes in FTS and coarse-
grained explicit-chain modelsPal2021 and have been applied more recently to
analyze experimental data on the biomolecular structure of nucleolar
condensates.JoshCliff The pair correlation functions are, in turn, related to
the potentials of mean force (PMFs), $U_{p,q}(r)$, through the standard
formula

$U_{p,q}(r)=-k_{\rm B}T\ln\left[\frac{G_{p,q}(r)}{\rho_{p}\rho_{q}}\right]\,,$
(51)

where the normalization factor $\rho_{p}\rho_{q}$ in the argument of the
logarithm is chosen such that $U_{p,q}<0$ means correlation and $U_{p,q}>0$
implies anti-correlation. In FTS, the pair correlation functions can be
computed using the following relations:

$\displaystyle\langle\hat{\rho}_{p}({\bf r})\hat{\rho}_{q}({\bf
r}^{\prime})\rangle$ $\displaystyle=\langle\tilde{\rho}_{p}({\bf
r})\tilde{\rho}_{q}({\bf r}^{\prime})\rangle\quad,\quad(p\neq q),$ (52)
$\displaystyle\langle\hat{\rho}_{p}({\bf r})\hat{\rho}_{p}({\bf
r}^{\prime})\rangle$ $\displaystyle=\left\langle\frac{\tilde{\rho}_{p}({\bf
r})}{v_{p}}\left(\frac{\mathrm{i}}{\gamma}\eta({\bf
r}^{\prime})+\rho_{0}-\sum_{q\neq p}v_{q}\tilde{\rho}_{q}({\bf
r}^{\prime})\right)\right\rangle\;.$

In accordance with a common procedure in field theory,IZ these relations can
be derived by adding source terms $\sum_{p}\int\mathrm{d}{\bf
r}\hat{\rho}_{p}({\bf r})J_{p}({\bf r})$ to the particle Hamiltonian $\hat{H}$
in Eq. (7) and then taking functional derivatives with respect to $J_{p}({\bf
r})$’s in the resulting field theory and finally setting all $J_{p}\rightarrow
0$. The formula for the self-correlation function
$\langle\hat{\rho}_{p}\hat{\rho}_{p}\rangle$ may be obtainedPal2021 through an
intermediate field re-definition to avoid a term that contains a double
derivative of $Q_{p}$.

Coarse-Grained Explicit-Chain Models. Following our previous studies,SumanPNAS
; Pal2021 ; wessen2021 coarse-grained explicit-chain molecular dynamics (MD)
simulations are conducted here to elucidate the ramifications of proposed
interaction schemesdignon18 ; SumanPNAS ; Urry-Mittal ; FB ; KH ; Mpipi and to
assess the accuracy of our approximate analytical theories and FTS in
capturing the essential physics of biomolecular LLPS. As before,dignon18 ;
suman2 each amino acid is represented by a single bead of different mass and
size. Although it is possible to incorporate explicit simple dipole solvent
molecules in coarse-grained explicit-chain LLPS simulation to better account
for dielectric effects,wessen2021 for simplicity we will limit the present
study to implicit-solvent simulations.dignon18 ; SumanPNAS As the general
features of the model can be found in the literature,dignon18 ; SumanPNAS ;
Urry-Mittal ; FB ; KH ; Mpipi only an outline is provided here. Briefly, as in
the above field-theory formulation, let $n_{\rm p}$ be the total number of
$N$-residue polypeptide chains (IDPs) in the MD system, with the chains
labeled by $i$ or $j=1,2,\dots,n_{\rm p}$ and the residues (beads) along a
polypeptide chain labeled by $\alpha$ or $\beta=1,2,\dots,N$. The total
potential energy $U_{\rm T}$ is a function of the residue positions $\\{{\bf
R}_{i,\alpha}\\}$ that may be written as

$U_{\rm T}=U_{\rm bond}+U_{\rm el}+U_{\rm hh}\;,$ (53)

where $U_{\rm bond}$ is the bond-length term for chain connectivity:

$U_{\rm bond}=\frac{K_{\rm bond}}{2}\sum_{i=1}^{n_{\rm
p}}\sum_{\alpha=1}^{N-1}(R_{i,\alpha+1;i,\alpha}-l)^{2}\;$ (54)

with $R_{i,\alpha;j,\beta}\equiv|{\bf R}_{i,\alpha}-{\bf R}_{j,\beta}|$.
Parameter values for the MD potential $U_{\rm bond}$ in Eq. (54) and for
$U_{\rm el}$, $U_{\rm hh}$ below are different across the three interaction
schemes---Mpipi,Mpipi Urry,Urry-Mittal and FBFB \---that we apply here to the
Ddx4 IDRs.SumanPNAS To faciliate comparisons, the parameters in these original
references are largely used with the corresponding interaction schemes in the
present simulations unless stated otherwise. For $U_{\rm bond}$, $K_{\rm
bond}=19.1917$ kcal mol-1Å-2, $l=3.81$ Å for Mpipi;Mpipi $K_{\rm bond}=9.56$
kcal mol-1Å-2, $l=3.82$ Å for Urry;Urry-Mittal and $K_{\rm bond}=2.39$ kcal
mol-1Å-2, $l=3.8$ Å for FB.FB The $U_{\rm el}$ term in Eq. (53) is the
electrostatic interaction:

$U_{\rm el}=\sum_{i,j=1}^{n_{\rm p}}\sum_{\begin{subarray}{c}\alpha,\beta=1\\\
\hbox{}\hskip-22.76228pt(i,\alpha)\neq(j,\beta)\end{subarray}}^{N}\frac{\sigma_{i,\alpha}\sigma_{j,\beta}e^{2}}{4\pi\epsilon_{0}\epsilon_{\rm
r}R_{i,\alpha;j,\beta}}\exp\Bigr{(}-\kappa_{\rm
D}R_{i,\alpha;j,\beta}\Bigr{)}\;,$ (55)

where, as in Eq. (13), $\sigma_{i,\alpha}$ is the electric charge, in units of
proton charge $e$, of the $\alpha$th residue along the $i$th chain
($\sigma_{i,\alpha}=\sigma_{\alpha}$ is independent of $i$ in the present MD
simulations), and $\kappa_{\rm D}$ is the reciprocal of the Debye screening
length. MpipiMpipi uses $\sigma_{\alpha}=+0.75$ for the positively charged
residues arginine and lysine, $\sigma_{\alpha}=+0.375$ for histidine (to
account for its incomplete ionization under physiological conditions), and
$\sigma_{\alpha}=-0.75$ for the negatively charged residues aspartic and
glutamic acids; UrryUrry-Mittal uses $\sigma_{\alpha}=+1$ for arginine and
lysine, $\sigma_{\alpha}=0$ for histidine, and $\sigma_{\alpha}=-1$ for
aspartic and glutamic acids; and FBFB uses $\sigma_{\alpha}=+1$ for arginine
and lysine, $\sigma_{\alpha}=+0.5$ for histidine, and $\sigma_{\alpha}=-1$ for
aspartic and glutamic acids. All other residues have $\sigma_{\alpha}=0$ in
all of these three interaction schemes. In all of the present coarse-grained
MD simulations, we set $\kappa_{\rm D}=0.1$ Å-1 (i.e., a Debye screening
length of 10 Å corresponding to a physiological monovalent salt concentration
of $\sim 100$ mM), apply a large distance cutoff of 35 Å to the electrostatic
interaction, and use a dielectric constant $\epsilon_{\rm r}=40$. Employing a
dielectric constant moderately lower than the $\epsilon_{\rm r}\approx 80$
value of bulk water in coarse-grained MD simulation of IDP LLPS has recently
been rationalized from a theoretical perspectiveSumanPNAS ; wessen2021 and
this modeling practice is apparently not inconsistent with experimental
observations of slower water orientational dynamics in living matter.Tros2017

For the spatially short-range interaction term $U_{\rm hh}$ in Eq. (53),
MpipiMpipi uses the Wang-Frenkel form,WangFrenkel viz.,

$U_{\rm hh}^{\rm Mpipi}=\sum_{i,j=1}^{n_{\rm
p}}\sum_{\begin{subarray}{c}\alpha,\beta=1\\\
\hbox{}\hskip-22.76228pt(i,\alpha)\neq(j,\beta)\end{subarray}}^{N}V_{\rm
WF}({\cal E}_{\alpha,\beta},\mu_{\alpha,\beta},1,(r_{\rm
c})_{\alpha,\beta},\sigma_{\alpha,\beta}|R_{i,\alpha;j,\beta})\,,$ (56)

where, by definition, $V_{\rm WF}({\cal
E}_{\alpha,\beta},\mu_{\alpha,\beta},1,(r_{\rm
c})_{\alpha,\beta},\sigma_{\alpha,\beta}$ vanishes for
$R_{i,\alpha;j,\beta}>(r_{\rm c})_{\alpha,\beta}$ as specified in Eqs. (47)
and (48). For all amino acid residue pairs $\alpha,\beta$, Mpipi sets the
distance cutoff $(r_{\rm c})_{\alpha,\beta}=3\sigma_{\alpha,\beta}$ and
assigns $\nu=1$. Values of the interaction strength ${\cal E}_{\alpha,\beta}$,
length scale $\sigma_{\alpha,\beta}$, and the parameter $\mu_{\alpha,\beta}$
for the Wang-Frenkel functional form depend on the amino acid residue types
$r_{\alpha}$ and $r_{\beta}$ for residues $\alpha$ and $\beta$. The Mpipi
values for ${\cal E}_{r,r^{\prime}}$, $\sigma_{r,r^{\prime}}$, and
$\mu_{r,r^{\prime}}$, where $r$ and $r^{\prime}$ each represents one of the
twenty amino acid types, are provided in Supplementary Table 11 of ref. (111),
wherein $\mu_{r,r^{\prime}}=2$ for almost all residue type pairs except for a
few cases in which $\mu_{r,r^{\prime}}=4$ or $11$.

For UrryUrry-Mittal and FBFB interaction schemes, $U_{\rm hh}$ in Eq. (53)
takes the following form:dignon18

$U_{\rm hh}^{\rm Urry/FB}=\sum_{i,j=1}^{n_{\rm
p}}\sum_{\begin{subarray}{c}\alpha,\beta=1\\\
\hbox{}\hskip-22.76228pt(i,\alpha)\neq(j,\beta)\end{subarray}}^{N}(U_{\rm
hh})_{i,\alpha;j,\beta}^{\rm Urry/FB}$ (57)

where

$(U_{\rm hh})_{i,\alpha;j,\beta}^{\rm Urry/FB}=\begin{cases}V_{\rm LJ}({\cal
E}_{\alpha,\beta},\sigma_{\alpha,\beta}|R_{i,\alpha;j,\beta})+(1-\lambda_{\alpha,\beta}^{\rm
hh}){\cal E}_{\alpha,\beta}&{\rm if\ }R_{i,\alpha;j,\beta}\leq
2^{1/6}\sigma_{\alpha,\beta}\\\ \lambda_{\alpha,\beta}^{\rm hh}V_{\rm
LJ}({\cal E}_{\alpha,\beta},\sigma_{\alpha,\beta}|R_{i,\alpha;j,\beta})&{\rm
otherwise}\end{cases}\;$ (58)

is in the Ashbaugh-Hatch formAH with the function $V_{\rm LJ}$ given by Eq.
(46). For both Urry and FB, ${\cal E}_{\alpha,\beta}={\cal E}=0.2$ kcal mol-1
irrespective of the residue types $r_{\alpha},r_{\beta}$ of residues
$\alpha,\beta$, and the length scale
$\sigma_{\alpha,\beta}=[\sigma(r_{\alpha})+\sigma(r_{\beta})]/2$, where
$\sigma(r_{\alpha})$ is the diameter of the bead model for amino acid residue
type $r_{\alpha}$ given by the quantity $\sigma$ in Table S1 of ref. (40). In
both interaction schemes, $\lambda_{\alpha,\beta}^{\rm
hh}=[\lambda(r_{\alpha})+\lambda(r_{\beta})]/2$, where $\lambda(r_{\alpha})$
is the hydrophobicity/hydropathy parameter for residue type $r_{\alpha}$. The
$\lambda(r_{\alpha})$ values for the Urry interaction scheme are given by the
‘‘Urry et al. normalized hydropathy scale’’ column of Table S2 in ref. (106),
whereas the corresponding $\lambda(r_{\alpha})$ values for the FB interaction
scheme are provided in Table S7 of ref. (107). In simulations that use the
Urry and FB interaction schemes, a cutoff distance of 20 Å is applied to the
$(U_{\rm hh})_{i,\alpha;j,\beta}^{\rm Urry/FB}$ interactions in Eq. (58).

In addition to modeling Ddx4 IDRs, to investigate nonelectrostatic effects of
sequence pattern on LLPS, coarse-grained explicit-chain MD simulations are
also conducted here for copolymer sequences consisting of only two types of
residues with zero electric charge but different hydrophobicities. The
hydrophobic-polar patterns of these model sequences---commonly referred to as
HP sequences in the protein literatureirback2020 ; dill1985 ; lau1989 ;
panag1992 \---are taken from recent studies of their LLPS.Statt2020 ;
panagio2021 In keeping with tradition of surfactant models, however,
hydrophobic and polar beads are labeled, respectively, as ‘‘T’’ (tail) and
‘‘H’’ (head) in refs. (126; 127). Here we use leucine for the hydrophobic
beads and serine for the polar beads. The resulting amino acid sequences are
thus referred to as LS sequences. As an example, we use the KH interaction
schemedignon18 described in ref. (66) for LS-sequence simulations. The
spatially short-range interaction $U_{\rm hh}$ is in the form of Eqs. (57) and
(58) wherein ${\cal E}_{\alpha,\beta}=0.228|(e_{\rm
MJ})_{r_{\alpha},r_{\beta}}-e_{0}|$, $\lambda^{\rm hh}_{\alpha,\beta}=1$ if
$(e_{\rm MJ})_{r_{\alpha},r_{\beta}}\leq e_{0}$ and $\lambda^{\rm
hh}_{\alpha,\beta}=-1$ otherwise, $e_{0}=-1.0$ kcal mol-1 and $(e_{\rm
MJ})_{r_{\alpha},r_{\beta}}$ is from ref. (109). The KH-D parameter set in
Table S3 of ref. (40) corresponds to $-\lambda^{\rm hh}_{\alpha,\beta}{\cal
E}_{r_{\alpha},r_{\beta}}$.

MD Simulation of LLPS. Based on the models described above, we follow the
recently developed, widely-applied slab method for simulation of IDP
LLPS.dignon18 ; suman2 ; mimb2022 Because this simulation protocol, its
rationale, and the procedure to construct phase diagrams from equilibrated
simulations have been detailed elsewhere,dignon18 ; suman2 ; panag2017 a brief
outline here will suffice. As in our previous studies,suman2 ; Pal2021 ;
SumanPNAS ; mimb2022 ; wessen2021 the GPU version of HOOMD-Blue softwareHOOMD
; Anderson is employed in the present simulations. For each simulation of the
wildtype or a variant of the 236-residue Ddx4 IDR,SumanPNAS 100 chains are
randomly placed in a sufficiently large cubic box of dimension $300\times
300\times 300$ Å3 initially. Energy minimization is then performed to remove
any steric clashes among the model molecules. This is followed by $NPT$
compression at a temperature of 100 K and pressure of 1 atm for a period of 50
ns using the Martyna-Tobias-Klein (MTK) thermostat and barostatklein1994 ;
martyna2006 with a coupling constant of 1 ps. The equations of motion are
integrated with velocity-Verlet algorithm using a large timestep of 10 fs.
Periodic boundary conditions are applied in all three directions. The
electrostatic interaction is treated with the PPPM algorithm.LeBard After the
initial $NPT$ compression, the system is compressed again at 100 K for a
period of 50 ns using Langevin thermostat with a friction coefficient of 1
ps-1 to reach a sufficiently high density phase of IDPs enclosing in a
simulation box of size $150\times 150\times 150$ Å3. At this point, the system
is expanded along one of the axes (referred to as the $z$-axis) to 3000 Å by
applying isotropic linear scaling for 10 ns with the temperature fixed at 100
K. Next, an $NVT$ equilibration is performed for 2 $\mu$s at desired
temperatures using the Langevin thermostat with a friction coefficient of 1
ps-1. Final production run is then carried out for another 4 $\mu$s with the
same Langevin thermostat using a lower friction coefficient of 0.01 ps-1.
Snapshots are saved every 1 ns for analysis. A similar procedure is used for
the simulations of the 20-residue LS (leucine-serine) sequences LS1, LS2, and
LS3. Here, the only difference with the Ddx4 IDR simulations is that 1,000 LS
chains are used for each simulation and the simulation box is compressed to a
size of $130\times 130\times 130$ Å3 and then expanded to 2000 Å along the
$z$-axis. Details regarding how to construct phase diagrams from equilibrium
trajectories are described in refs. (40; 43; 66; 78).

RESULTS AND DISCUSSION

While we consider the above-detailed theoretical development to be the main
thrust of this article and that extensive testing of our theory is beyond the
scope of the present work, it is instructive to apply the new formulation to a
few initial examples to probe our theory’s practical effectiveness and to
provide suggestions for future theoretical/computational improvements from the
success and limitation of these examinations. We do so below by using our
formulation on four well-studied Ddx4 IDRsNott15 ; jacob2017 ; SumanPNAS ;
robert ; Urry-Mittal ; Mpipi as well as three selected hydrophobic-polar model
sequences that have the same composition but different hydrophobic/polar
sequence patterns.Statt2020

Analytical Theory for the Phase Behaviors of Wildtype and Variant Ddx4 IDRs
Modeled by Different Interaction Schemes. To illustrate the approximate
analytical approach in the previous section, we apply the theory with the free
energy function given by Eq. (34) to the 236-residue wildtype (WT), charge-
scrambled (CS), arginine-to-lysine (RtoK), and phenylalanine-to-alanine (FtoA)
variants of the Ddx4 IDR (the amino acid sequences are provided, for example,
in Fig. S1 of ref. (66)). We compare the LLPS properties predicted by using
$\varepsilon_{r,r^{\prime}}$s from the five interaction schemes illustrated in
Fig. 2. In accordance with Eqs. (36) and (37) for $\chi_{\rm eff}$, all
nonzero eigenvalues $\lambda_{a}$s are included in the application of the
MFT/RPA Eq. (34) to the four Ddx4 IDRs.

Fig. 5: Nonelectric charge profiles for Ddx4 IDRs. Distributions of
nonelectric charges $q^{(a)}{}_{r_{\alpha}}$ for the two most dominant
eigenvalues $a=1,2$ of $\varepsilon_{r,r^{\prime}}$ [Eq. (1)] along the WT,
RtoK and FtoA Ddx4 sequences are shown in the (a) KH, (b) Mpipi, (c) HPS, (d)
Urry, and (e) FB interaction schemes [key for line styles provided in (a)].
The nonelectric charge profiles in this figure are averages of individual
residues’ $q^{(1)}{}_{r_{\alpha}}$ or $q^{(2)}{}_{r_{\alpha}}$ over a sliding
window of six residues. As generally defined in the text,
$\alpha=1,2,\dots,236$ is the sequential position label of the residues along
the Ddx4 IDRs, and $r_{\alpha}$ is the amino acid type of residue $\alpha$.

Fig. 5 shows the nonelectric charge profiles of the WT, RtoK and FtoA Ddx4
IDRs for the two most dominant eigenvalues $\lambda_{1,2}$ in each of the five
interaction schemes considered. Significant differences in nonelectric charge
profiles are observed across different interaction schemes (Fig. 5a--e),
reflecting variations in the assumptions made by KH, Mpipi, HPS, Urry, and FB
regarding physical driving forces in biomolecular processes. Variation in the
two dominant nonelectric charges along the IDR sequences (as function of the
horizontal varible $\alpha$) is discernibly more pronounced in Mpipi and FB
than in HPS and Urry, with KH exhibiting an intermediate degree for variation
along the chain sequence. Implications of this difference for predicted LLPS
propensity remain to be ascertained. Within the same interaction scheme, the
variation in nonelectric charge profiles among the three Ddx4 sequence
variants (shown by curves in different colors in Fig. 5) is notably higher for
the KH and Mpipi interaction schemes with 210 contact energies than for the
HPS, Urry, and FB interaction schemes that are based upon 20-value
hydrophobicity/hydropathy scales, indicating that the LLPS propensities of the
sequence variants predicted by KH and Mpipi are likely more diverse in KH and
Mpipi than those predicted by HPS, Urry, and FB.

Because the WT, RtoK and FtoA Ddx4 IDRs share the same sequence electric
charge pattern, any sequence dependence of LLPS propensity among these three
sequences in the MFT/RPA analytical approach originates from the difference in
their $\chi_{\rm eff}$ parameters, which may be computed from the eigenvalues
$\lambda_{a}$ of the spectral decomposition of $\varepsilon_{r,r^{\prime}}$ in
conjunction with the nonelectric charges $q^{(a)}{}_{r_{\alpha}}$ using Eq.
(36) or directly from the interaction matrix $\varepsilon_{r,r^{\prime}}$
itself via the identity in Eq. (37). In the present MFT/RPA formulation, which
for simplicity assigns the same $V_{\rm h}(r)$ to every residue-residue
interaction scheme considered, the ratios between these $\chi_{\rm eff}$
parameters depend only on $\varepsilon_{r,r^{\prime}}$ but not $V_{\rm h}(r)$.
As such, these ratios may serve as a zeroth order approximate measure of the
relative contributions to LLPS propensity from the spatially short-range,
nonelectrostatic interactions encoded by the different IDR sequences.
Accordingly, for the five interaction schemes we consider, the relative short-
spatial-range interaction strengths among the WT, CS, RtoK and FtoA Ddx4 IDRs
is quantified by MFT in terms of the ratios

$\chi^{\rm(WT,CS)}_{\rm eff}:\chi^{\rm(RtoK)}_{\rm eff}:\chi^{\rm(FtoA)}_{\rm
eff}\;$ (59)

wherein WT and CS share the same $\chi_{\rm eff}$ because the WT and CS
sequences have the same amino acid composition. Utilizing the
$\varepsilon_{r,r^{\prime}}$ values of the interaction schemesdignon18 ; Urry-
Mittal ; FB ; Mpipi in Eqs. (36) and (37) leads to the following ratio:

$\displaystyle\chi^{\rm(WT,CS)}_{\rm eff}:\chi^{\rm(RtoK)}_{\rm
eff}:\;\chi^{\rm(FtoA)}_{\rm eff}$ (60) $\displaystyle\quad\quad\quad
1\,:\;\;0.83\;\,:\;\,0.66\quad$ $\displaystyle(\mbox{KH})$
$\displaystyle\quad\quad\quad 1\,:\;\;0.80\;\,:\;\,0.84\quad$
$\displaystyle(\mbox{Mpipi})$ $\displaystyle\quad\quad\quad
1\,:\;\;1.09\;\,:\;\,0.97\quad$ $\displaystyle(\mbox{HPS})$
$\displaystyle\quad\quad\quad 1\,:\;\;0.96\;\,:\;\,0.97\quad$
$\displaystyle(\mbox{Urry})$ $\displaystyle\quad\quad\quad
1\,:\;\;1.04\;\,:\;\,0.94\quad$ $\displaystyle(\mbox{FB})\;.$

Fig. 6: Salt-dependent phase behaviors of Ddx4 IDRs computed by analytical
theory. Phase diagrams of the four Ddx4 variants (color-coded on the right) on
the protein-salt density $(\rho_{\rm p},\rho_{\rm s})$ plane are computed in
our approximate MFT/RPA analytical approach [Eq. (34)] at model temperature
$T=280$ K using five different interaction schemes as indicated by the top
labels. The phase-separated regimes are enclosed by coexistence (binodal)
curves (solid curves), $(\rho_{\rm p},\rho_{\rm s})$ of coexisting phases are
connected by tie-lines (dashed lines with the same color code for the sequence
variants). All parameter values used for this figure are described in the
text.

Phase diagrams can readily be constructedmimb2022 using the MFT/RPA free
energy function in Eq. (34). To facilitate comparison of the five different
interaction schemes, the overall strength of the spatially short-range
interactions, encapsulated by $\hat{V}_{\rm h}(\bm{0})$, is adjusted for each
interaction scheme separately such that $\chi^{\rm(WT,CS)}_{\rm eff}=0.5$ at a
reference temperature of $T=300$ K. We also set $l_{\rm B}=7$ Å at the same
reference temperature throughout, and assume the temperature dependence of
$l_{\rm B},\hat{V}_{\rm h}(\bm{0})\propto T^{-1}$ (as stated above, $V_{\rm
h}(r)$ is potential energy in units of $k_{\rm B}T$). Other parameters used
are $v_{\rm w}=1$, $v_{\rm b}=2$, $v_{\pm}=0$, $a_{\rm s}=b/\sqrt{6}$, and
$\rho_{0}=55.5$ M.

Fig. 6 shows phase diagrams in the $(\rho_{\rm p},\rho_{\rm s})$ plane, where
the two-phase region is enclosed by the coexisting (binodal) curves (solid
lines). Systems with bulk concentrations within the two-phase region separate
into two phases with the partition of molecular species provided by the tie-
lines (dashed lines in Fig. 6). The phase diagrams in Fig. 6 are computed at a
uniform temperature of $T=280$ K, which is chosen to be sufficiently low such
that every Ddx4 IDR phase separates under any one of the spatially short-range
interaction scheme we consider.

Fig. 7: Temperature-dependent phase behaviors of Ddx4 IDRs computed by
analytical theory. As in Fig. 6, phase diagrams are computed in our
approximate MFT/RPA analytical approach [Eq. (34)] using five different
interaction schemes. Shown here are protein density ($\rho_{\rm p}$) versus
model temperature ($T$) phase diagrams for overall (bulk) monovalent salt
([NaCl]) concentrations of 100 mM (top row) and 300 mM (bottom row). These are
constructed using constant-$T$ $(\rho_{\rm p},\rho_{\rm s})$ phase diagrams
(exemplified by Fig. 6) at multiple $T$. For a given $T$, the two coexisting
$\rho_{\rm p}$s at bulk salt density $\rho_{\rm s}=100$ or $300$ mM are
obtained as the $\rho_{\rm p}$ values at the two ends of a $(\rho_{\rm
p},\rho_{\rm s})$ coexistence tie-line that passes through the intersection
point of $\rho_{\rm p}=10$ mM (indicated by the vertical dotted lines in the
present figure) and $\rho_{\rm s}=100$ mM or $\rho_{\rm s}=300$ mM,
respectively, in the $\rho_{\rm p}$\--$\rho_{\rm s}$ plane. The color code for
the Ddx4 variants are the same as that in Fig. 6. Dashed curves are use here
for FtoA to facilitate visual discrimination between the near-overlapping
coexisting curves for FtoA (red) and RtoK (green) in the Urry interaction
scheme.

Fig. 7 shows $(\rho_{\rm p},T)$ phase diagrams of the same systems at a
constant bulk protein density of $\rho_{\rm p}=10$ mM (indicated by the
vertical grey dashed lines), and at two bulk salt densities $\rho_{\rm s}=100$
mM (top) and 300 mM (bottom). The phase diagrams were constructed by first
computing constant-temperature binodal curves (such as Fig. 6) in the range
$150\mbox{ K}\leq T\leq 510\mbox{ K}$. At each temperature, we check if the
point $(\rho_{\rm p},\rho_{\rm s})$ lies in the two-phase region and then
identify the tie-line it intersects. The densities of the co-existing phases
are then given by the two points on the binodal curve connected by the tie-
line. It should be noted however that, by construction, the coexisting protein
densities provided by Fig. 7 for a given temperature and a given bulk salt
concentration apply only to systems with bulk protein density $\rho_{\rm
p}=10$ mM. The coexisting protein densities (i.e., the phase diagram) for the
same temperature and bulk salt concentration would be different for systems
with a different bulk protein density. In this regard, the meaning of the
phase diagrams in Fig. 7 is more restrictive than phase diagrams for one
solute species in which the coexisting solute densities for a given
temperature apply to all bulk solute densities within the two-phase region.

Results in Figs. 6 and 7 confirm and add to previous observations that LLPS
properties entailed by different interaction schemes can be significantly
different.SumanPNAS ; Urry-Mittal ; Mpipi Two-solute $(\rho_{\rm p},\rho_{\rm
s})$ phase properties of Ddx4 IDRs such as those depicted in Fig. 6 have not
been considered before. A novel feature that emerges in Fig. 6 is that, for
the KH and Mpipi interaction schemes, the coexistence curve for the CS variant
(orange curve) intersects the coexistence curves of other Ddx4 variants (green
and red curves), indicating that the rank ordering of LLPS propensities is
predicted by these interaction schemes to be salt dependent. This is caused by
the fact that the sequence electric charge pattern of CS is different from the
other Ddx4 variants and therefore the electrostatic screening effects of salt
on CS phase separation vary in a different manner from that on the other Ddx4
variants when salt concentration varies. The same feature is manifested by the
crossing between the CS coexistence curves and those of other variant(s) in
the KH and Mpipi panels in Fig. 7, indicating that the rank ordering of LLPS
propensities of the Ddx4 IDRs is predicted by these two interaction schemes to
be temperature dependent as well.

As noted before,SumanPNAS ; Mpipi for KH and Mpipi, all three variants clearly
exhibit lower LLPS propensities than the WT (the blue curve is farther to the
right than the other curves in the KH and Mpipi panels in Fig. 6, and the blue
curve for the condensed-phase protein density is higher than the other curves
in the corresponding panels in Fig. 7). This trend is consistent with
experiments.Nott15 ; jacob2017 ; robert For the RtoK and FtoA variants, the
trend is underpinned by the effective Flory $\chi$ parameter ratios
$\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.83$,
$\chi^{\rm(FtoA)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.66$ for KH and
$\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.80$,
$\chi^{\rm(FtoA)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.84$ for Mpipi in Eq.
(60). In the Urry interaction scheme, the three Ddx4 variants also show lower
LLPS propensities than that of the WT, but only barely for RtoK and FtoA
($\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.96$,
$\chi^{\rm(FtoA)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.97$ for Urry). In
contrast, as pointed out before,SumanPNAS results for HPS in Figs. 6 and 7
indicate that the RtoK variant has a higher LLPS propensity than WT, which is
inconsistent with experimentsNott15 ; jacob2017 ; robert
($\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=1.09$,
$\chi^{\rm(FtoA)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.97$ for HPS). The
results in Figs. 6 and 7 for FB also exhibit a higher LLPS propensity for RtoK
than for WT, though to a lesser degree. Similar to HPS, this FB prediction
($\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=1.04$,
$\chi^{\rm(FtoA)}_{\rm eff}/\chi^{\rm(WT)}_{\rm eff}=0.94$ for FB) is caused
by its exceptionally high hydrophobicity parameter of $0.47106$ for lysine
compared to that of $0.24025$ for arginine (Table S7 for optimized FB-HPS
parameters in ref. (107)).

Experimental data on Ddx4 IDRs indicate that LLPS propensity of WT is higher
than that of CS (ref. (9)) and that LLPS propensity of CS is higher than that
of FtoA (ref. (19), FtoA corresponds to 14FtoA in this reference). Because
LLPS for FtoA is observed at $\sim 350$ mg mL-1 protein concentration but not
for RtoK up to $400$ mg mL-1 protein concentration under the conditions given
in Appendix 1--Table 4 of ref. (96)), the overall experimental rank ordering
of LLPS propensities is WT $>$ CS $>$ FtoA $>$ RtoK. Notably, this rank
ordering is reproduced by Mpipi, as shown in Fig. 6b of ref. (111) and Figs. 6
and 7 here (red curves are higher than green curves for Mpipi), with
$\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(FtoA)}_{\rm eff}=0.95$ for Mpipi. By
comparison, KH does not capture this trend entirely, as its predicted rank
ordering of LLPS propensities is WT $>$ CS $>$ RtoK $>$ FtoA
($\chi^{\rm(RtoK)}_{\rm eff}/\chi^{\rm(FtoA)}_{\rm eff}=1.26$ for KH), as
shown in Fig. 4 of ref. (66) and in Figs. 6 and 7 here (red curves are lower
than green curves for KH).

Field-Theoretic Simulations of Phase Behaviors of Model Hydrophobic-Polar
Sequences with the Same Composition but Different Sequence Patterns. As a test
of the FTS approach introduced above for modeling sequence-dependent LLPS
effects of spatially short-range interaction, we apply the formulation to
three model hydrophobic-polar sequences with the same $3/2$ hydrophobic/polar
composition but different sequence arrangements of the hydrophobic and polar
residues. We adopt three 20mer ($N=20$) sequences with $f_{\rm T}=0.6$ from
Statt et al.Statt2020 (hydrophobic composition is denoted as $f_{\rm T}$ in
this reference), namely their ‘‘[HT]2TH[T2H]4TH’’, ‘‘T3H3T3H2T3H3T3’’, and
‘‘T12H8’’ sequences listed in ref. (126) as, respectively, the last sequence
in Fig. S6, the sequence in Fig. 5, and the ‘‘micelles’’ sequence in Fig. 6.
As mentioned above, since we use leucine and serine, respectively, for the
hydrophobic and polar residues in the present study, we refer to these three
sequences as LS1, LS2, and LS3 (Fig. 8a). The hydrophobic and polar residues
are more evenly distributed in LS1 and LS2 than in the diblock arrangement of
LS3. Statt et al. essentially use a standard hydrophobic-polar (HP)
potentiallau1989 whereby only hydrophobic-hydrophobic interactions are
favorable (governed by a LJ potential [Eq. (46) with negative contact energy
and excluded volume], polar-polar and polar-hydrophobic interactions are
repulsive because of excluded volume (governed by a WCA potentialWCA ). Within
this HP interaction scheme, the LLPS critical temperatures ($T_{\rm cr}=T_{c}$
in their notation) for the first and second sequences are, in model units,
$T_{\rm cr}=0.808$ and $T_{\rm cr}=1.083$, whereas the third sequence
apparently does not have a critical point (Table S1 of ref. (126)). Since we
will be applying the KH interaction scheme to the LS sequences, with ${\cal
E}_{r_{\alpha},r_{\beta}}=-0.771,0$, and $-0.304$, respectively, for
$(r_{\alpha},r_{\beta})=$ (leucine, leucine), (serine, serine), (leucine,
serine), instead of the HP potential, LLPS behaviors of our LS sequences are
expected to exhibit a similar but not identical trend to their corresponding
HP-like sequences.

Fig. 8: Sequence-specific LLPS effects of hydrophobic-like short-spatial-range
interactions. (a) Leusine-Serine (LS) 20mer polymer sequences used for FTS in
this work. (b,c) Polymer-polymer and polymer-solvent pair correlation
functions $G_{\rm b,b}(r)$ and $G_{\rm b,s}$, defined in Eq. (50), computed in
FTS ($p={\rm b}$, $q={\rm s}$). (d,e) PMFs computed from the pair-correlation
functions in (a,b) using Eq. (51). Solid curves and shaded bands in panels (b
--e) depict the average and standard deviation of 21 independent simulations.
The gray solid curve in each panel corresponds to $l_{\rm h}=0.03$ Å, wherein
short-spatial-range interactions are not sufficiently strong to drive
condensate formation and for which no visible sequence dependence is observed.
Coloured curves are computed at $l_{\rm h}=0.5$ Å for which the strong short-
spatial-range interactions cause the polymers to form a dense condensate.

Our FTS simulations are performed on an $M=48$ lattice with lattice spacing
$\Delta x=a_{\rm s}=b/\sqrt{6}\approx 1.6$ Å, which amounts to a simulation
box width of around $75$ Å. Using the KH interaction matrix, we truncate the
summations over eigenvalues at $\tilde{a}=3$ [Eqs. (1) and (2)]. The range of
spatially short-range interactions governed by $V_{\rm h}(r)$ in Eq. (49) is
set to $\kappa^{-1}=5$ Å and we let $v_{\rm b}=v_{\rm w}=1$ for simplicity.
All simulations are performed without salt, and since the LS sequences carry
no electric charge, the field $\psi({\bf r})$ in Eq. (20) can be omitted. The
total density of the system is set to $\rho_{0}=\rho_{\rm w}+N\rho_{\rm
p}=55.5$ M, with a polymeric (LS sequence) component of $\rho_{\rm p}=300$ mM
giving a polymer bead volume fraction of $\phi=0.11$. The system thus contains
$n_{\rm p}\approx 75$ LS polymer chains and $n_{\rm w}\approx 1.2\cdot 10^{4}$
solvent molecules. We set the compressibility to $\gamma=0.1b^{-3}\approx 3.0$
M because of equilibration issues observed in FTS of exactly incompressible
systems.wessen2021 However, we note that incompressible systems can
alternatively be explored in FTS using the partial saddle-point approximation
of ref. (155). Simulations are performed at $l_{\rm h}=0.5$ Å ($\chi_{\rm
eff}=2.16$) and at $l_{\rm h}=0.03$ Å ($\chi_{\rm eff}=0.129$), corresponding,
respectively, to strong and weak short-spatial-range interactions. The CL
evolution equations in Eq. (38) are numerically solved with a time-step
$\Delta t=5\cdot 10^{-4}b^{3}$ using a semi-implicit integration
schemeLennon2008 generalized to multiple fields.mimb2022 We perform 21
independent simulations for each sequence and parameter set, and use $Q_{\rm
p,w}$ to monitor the system equilibration. At $l_{\rm h}=0.5$ Å, the system
requires $4\cdot 10^{5}$ timesteps for equilibration, after which the field
configuration are sampled every 50th step during the production run consisting
of another $4\times 10^{5}$ CL steps. At $l_{\rm h}=0.03$ Å, $2\times 10^{4}$
integration steps are sufficient for equilibration, followed by sampling every
50th step during a production run of $3\times 10^{5}$ steps.

Fig. 8b--e shows pair-correlation functions and PMFs computed in FTS using
Eqs. (50) and (51) for the three LS sequences. The grey curve in each panel
corresponds to $l_{\rm h}=0.03$ Å, for which the system remains in a
homogeneous (single-phase) state exhibiting no significant sequence dependence
in the properties simulated. This single-phase state is characterized by all
PMFs approaching zero at large $r$. Coloured curves in Fig. 8 correspond to
$l_{\rm h}=0.5$ Å, at which the short-spatial-range interactions are
sufficiently strong to drive the system into an inhomogeneous state containing
a single spherical polymer-dense condensate against a near polymer-empty
background. This feature is characterized by the $G_{\rm p,p}(r)\sim 0$ (Fig.
8b) at separations $r$ larger than the condensate diameter, combined with
polymer-solvent anti-correlation at small $r$ because solvent are being pushed
out of the condensate.

Sequence-dependent effects are exhibited by the correlation functions and the
PMFs. While the sequences LS1 and LS2 behave very similarly (blue and orange
curves essentially overlap in Fig. 8b--e), the diblock sequence LS3 stands out
clearly by forming denser condensates than both LS1 and LS2 (at small $r$, the
red curve is higher than the orange and blue curves for correlation functions
in Fig. 8b; and the red curve is lower than the orange and blue curves for
PMFs in Fig. 8d). We also find that the LS3 condensate possesses a more
heterogeneous internal structure with distinct sub-regions populated
predominantly by either L or S residues. This feature can be seen in the
simulation snapshots in Figs. 9 and 10, which depict the real component of the
field operators

$\tilde{\rho}_{X}({\bf r})=\frac{\rho_{\rm p}}{Q_{\rm
p}}\sum_{\alpha=1}^{N}\mathrm{e}^{W_{\alpha}({\bf r})}q_{\rm F}({\bf
r},\alpha)q_{\rm B}({\bf r},\alpha)\delta_{r_{\alpha},X},\quad X=\mbox{L
(leusine)},\mbox{S (serine)},$ (61)

for the number density of L and S residues (the Kronecker
$\delta_{r,r^{\prime}}=1$ for $r=r^{\prime}$ and $\delta_{r,r^{\prime}}=0$ for
$r\neq r^{\prime}$). To facilitate visualization of the condensate structure,
the snapshots in Fig. 9 are averaged over $10^{3}$ consecutive CL time steps.
Fig. 10 provide cross-sectional views through the center of mass of each of
the snapshots in Fig. 9. Leucine and serine densities are represented by
contours in Fig. 10 following the style introduced by Fig. 3 of ref. (62). The
contour plots in Fig. 10 demonstrate a marked increased degree of structural
heterogeneity of the LS3 condensate relative to the LS1 and LS2 condensates.
For LS3, the serine (polar) residues are predominantly on the outside whereas
the leucine (hydrophobic) residues are predominantly inside. This phenomenon
is reminiscent of micellar structure and, in this regard, consistent with the
simulated structure of the corresponding hydrophobic-polar sequence in ref.
(126).

Fig. 9: Field snapshots illustrate LLPS of LS sequences. Snapshots of FTS
after equilibration, showing field operators corresponding to leusine and
serine number densities in orange and cyan, respectively. The left, middle and
right columns are for the sequences LS1, LS2 and LS3, respectively, in Fig.
8a. (a) Top panels show snapshots under conditions where interactions are not
sufficiently strong to cause condensation ($T>T_{\rm cr}$, where $T_{\rm cr}$
is upper critical solution temperature, UCST). (b) Bottom panels depict
systems under conditions where there are sufficiently strong interactions for
phase separation ($T<T_{\rm cr}$). All snapshots are averaged over $10^{3}$ CL
steps.  Fig. 10: Internal structure of condensates of model hydrophobic-polar
sequences. Cross-sectional density distributions of the snapshots in Fig. 9 in
planes indicated by the insets. Leusine and serine densities are shown in
orange and cyan, respectively, as in Fig. 9. (a) Snapshots at $l_{\rm h}=0.03$
Å where the systems remain in a single-phase state. (b--d) Snapshots at
$l_{\rm h}=0.5$ Å where the LS polymers coalesce into a near-spherical
condensate.

Assessing Analytical and FTS Results by Coarse-Grained Explicit-Chain MD. We
now proceed to further evaluate the effectiveness of above analytical and FTS
approaches by comparing their predictions with the corresponding results from
coarse-grained explicit-chain models. As described in Models and Methods,
while field-theoretic formulation and coarse-grained modeling aim to capture
essential physics of the same system, details of the two types of models
differ. For instance, functional form of residue-residue interactions is
practically limited in field-theoretic formulation because of mathematical
consideration---in the present case we are restricted to a single Yukawa form
for all residue pairs, i.e., using $V_{\rm h}(r)$ in Eq. (49) with the same
$l_{\rm h}$ and $\kappa$ for all amino acid pairs. In contrast, there is more
freedom in selecting functional forms for potential energy (e.g., Lennard-
Jones, Ashbaugh-Hatch, Wang-Frenkel) as well as using different bead sizes to
represent different amino acid residues; e.g., different
$\sigma_{\alpha,\beta}$ for different types of amino acid pairs
$r_{\alpha},r_{\beta}$ in Eq. (56). This recognition notwithstanding, to be
useful as a computationally efficient approach complementary to coarse-grained
explicit-chain MD, field-theoretic method applied to any given system should
produce results that are in agreement, at least semi-quantitatively, with
those obtained from an appropriately constructed coarse-grained explicit-chain
model for the same system.

Coarse-Grained Explicit-Chain MD for Ddx4 IDRs. With this in mind, phase
diagrams for the four Ddx4 IDRs are obtained by coarse-grained explicit-chain
MD using the Mpipi,Mpipi Urry,Urry-Mittal and FBFB interaction schemes (Fig.
11). Corresponding coarse-grained explicit-chain MD phase diagrams using KH
and HPS have been provided, respectively, in Figs. 4 and 3B of ref. (66). The
trends exhibited by these coarse-grained MD results are largely consistent
with the analytical MFT/RPA results in Figs. 6 and 7, suggesting that the
MFT/RPA formulation is a computationally efficient approach for exploring
sequence-dependent LLPS of IDPs. In particular, in agreement with Fig. 6b of
the original Mpipi studyMpipi (which uses $\kappa_{\rm D}^{-1}=7.95$ Å,
$\epsilon_{\rm r}=80$ instead of $\kappa_{\rm D}^{-1}=10$ Å, $\epsilon_{\rm
r}=40$ in Fig. 11a here) and consistent with the present analytical MFT/RPA
results in Fig. 7 at lower model temperatures $T\lesssim 250$ K, Fig. 11a
exhibits the experimentally correct WT $>$ CS $>$ FtoA $>$ RtoK rank ordering
of LLPS propensities,Nott15 ; jacob2017 ; robert though the Mpipi-computed
difference in LLPS propensity between WT and CS ($\sim 6$\--$10$ K difference
in $T_{\rm cr}$)Mpipi is significantly smaller than the experimental
estimation of $\sim 80$ K difference in $T_{\rm cr}$ at [NaCl] = 100 mM (ref.
(19)). For Urry, the LLPS propensity rank ordering WT $>$ FtoA $\approx$ RtoK
$>$ CS computed by coarse-grained MD (Fig. 11b) is consistent with that for
Urry in MFT/RPA Figs. 6 and 7. Similarly, for FB, the RtoK $>$ WT $>$ FtoA $>$
CS rank ordering in Fig. 11c is consistent with that for FB in Figs. 6 and 7.

Fig. 11: Comparing explicit-chain model phase behaviors of Ddx4 IDRs in
different interaction schemes. Phase diagrams (coexistence curves) for WT
(blue), CS (orange), FtoA (red), and RtoK (green) are obtained by coarse-
grained MD at $\epsilon_{\rm r}=40$ using the (a) Mpipi,Mpipi (b) Urry,Urry-
Mittal or (c) FBFB interaction schemes in accordance with the modeling details
and simulation procedure described in Models and Methods. Lines passing
through simulation data points and extrapolated critical pointsdignon18 ;
SumanPNAS ; panag2017 (filled circles) are a guide for the eye.

As discussed in previous works,SumanPNAS ; Urry-Mittal ; FB ; Mpipi some of
the salient differences among interaction schemes can be understood semi-
quantitatively by simple considerations of the schemes’ contact energies. For
instance, the differences in LLPS behavior of the RtoK variant across
different interaction schemes are clearly related to the average strength of
interactions involving R versus that involving K in the interaction schemes.
For Mpipi, consider the average ${\cal E}_{\rm R,K}\equiv\sum_{r=1}^{20}{\cal
E}_{r,r^{\prime}}/20$ over the 20 amino acid types, where $r^{\prime}$ is the
label for arginine (R) or lysine (K). Entries for ${\cal E}_{r,r^{\prime}}$ in
Supplementary Table 11 of ref. (111) yield ${\cal E}_{\rm R}=0.2212$ kcal
mol-1 and ${\cal E}_{\rm K}=0.0563$ kcal mol-1, thus $\Delta{\cal E}_{\rm
R-K}\equiv{\cal E}_{\rm R}-{\cal E}_{\rm K}=$ $0.1649$ kcal mol-1 may be used
to characterize the degree to which interactions involving R are more
favorable than those involving K. The corresponding quantity in KH is
$\Delta(\lambda^{\rm hh}{\cal E})_{\rm R-K}=0.0829$ kcal mol-1, which is
significantly smaller. Here, for KH, $\Delta(\lambda^{\rm hh}{\cal E})_{\rm
R-K}\equiv(\lambda^{\rm hh}{\cal E})_{\rm R}-(\lambda^{\rm hh}{\cal E})_{\rm
K}$, with $\lambda^{\rm hh}{\cal E}_{\rm
R,K}\equiv-\sum_{r=1}^{20}\lambda^{\rm hh}_{r,r^{\prime}}{\cal
E}_{r,r^{\prime}}/20$, and $-\lambda^{\rm hh}_{r,r^{\prime}}{\cal
E}_{r,r^{\prime}}$ being the entries in Table S3 of ref. (40). For HPS, Urry,
and FB, the corresponding difference in interaction strength is given by
${\cal E}\Delta\lambda_{\rm R-K}/2\equiv{\cal E}[\lambda({\rm R})-\lambda({\rm
K})]/2$ where ${\cal E}=0.2$ kcal mol-1 and $\lambda({\rm R})$ and
$\lambda({\rm K})$ are, respectively, the value for R and K on a given
hydrophobicity/hydropathy scale (see Models and Methods). From the scales in
Table S1 of ref. (40), the ‘‘Urry et al.’’ column of Table S2 of ref. (106),
and Table S7 of ref. (107), ${\cal E}\Delta\lambda_{\rm R-K}/2=$ $-0.0514$,
$0.0176$, and $-0.0231$ kcal mol-1, respectively, for HPS, Urry, and FB. Taken
together, these numbers for Mpipi, KH, Urry, HPS, and FB, viz., $0.1807$,
$0.0829$, $0.0176$, $-0.0459$, and $-0.0170$ respectively, provide a useful
explanation for the signs and magnitudes of their WT $-$ RtoK critical
temperature differences $\Delta(T_{\rm cr})_{\rm WT-RtoK}\equiv(T_{\rm
cr})_{\rm WT}-(T_{\rm cr})_{\rm RtoK}$ $\approx$ $60$, $43$, $4$, $-29$, and
$-24$ K as seen for Mpipi, KH, Urry, HPS, and FB, respectively, in Fig. 11 of
the present study as well as Figs. 3B and 4 of ref. (66).

As noted above, the difference in LLPS propensity between WT and CS predicted
by Mpipi is small, with $\Delta(T_{\rm cr})_{\rm WT-CS}\approx 6$ K in the
original simulationMpipi and $\approx 15$ K in Fig. 11a, both much less than
the experimentally estimated $\Delta(T_{\rm cr})_{\rm WT-CS}\sim 80$ K (refs.
(19)). In fact, the slightly higher $\Delta(T_{\rm cr})_{\rm WT-CS}$ we obtain
is probably due in large part to the smaller dielectric constant
$\epsilon_{\rm r}=40$ used for Fig. 11a instead of the $\epsilon_{\rm r}=80$
value used in the original study. By comparison, larger $\Delta(T_{\rm
cr})_{\rm WT-CS}$ values are predicted by most of the other interaction
schemes considered here: For HPS, $\Delta(T_{\rm cr})_{\rm WT-CS}\approx 14$ K
for $\epsilon_{\rm r}=80$ (Fig. 3B of ref. (66)). For KH, $\Delta(T_{\rm
cr})_{\rm WT-CS}\approx$ 21 and 46 K, respectively, for $\epsilon_{\rm r}=80$
and $40$ (Fig. 4 of ref. (66)). For Urry and FB simulated here at
$\epsilon_{\rm r}=40$, $\Delta(T_{\rm cr})_{\rm WT-CS}\approx 60$ K (Fig.
11b,c). Presumably, this conspicuous difference between Mpipi and the other
interaction schemes regarding CS properties is chiefly caused by Mpipi’s
electric charge assignment of $0.75e$ for R and K, $-0.75e$ for aspartic acid
(D) and glutamic acid (E), and $0.375e$ for histidineMpipi instead of using
the full proton charge $e$ for electrostatic interactions among R, K, D, and E
as in the other interaction schemes. Since each of the Ddx4 IDRs contains only
two histidines, the main effects are in R, K, D, and E. The assignment of $\pm
0.75e$ electric charges instead of $\pm e$ amounts to a reduction factor
$(0.75)^{2}=0.5625$ in electrostatic interaction strength, which is further
attenuated by Mpipi’s adoption of a slightly shorter Debye screening length
$\kappa_{\rm D}^{-1}=7.95$ Å instead of the $\kappa_{\rm D}^{-1}=10$ Å
employed for the other interaction schemes. Because the difference in LLPS
propensity between WT and CS is principally a sequence-specific electrostatic
effect, Mpipi’s reduced electrostatic strength should go a long way in
accounting for the relatively small $\Delta(T_{\rm cr})_{\rm WT-CS}$ it
predicts.

Secondary to electrostatics, another possible origin of the observed variation
in $\Delta(T_{\rm cr})_{\rm WT-CS}$ predicted by different interaction schemes
is their different spatially short-range nonelectrostatic residue-residue
contact energies. WT is transformed into CS, and vice versa, by swapping
residues among 24 sequence positions---involving 6 D, 6 E, 9 R, and 3 K
residues participating in 4 (E $\leftrightarrow$ R), 5 (D $\leftrightarrow$
R), 2 (E $\leftrightarrow$ K), and 1 (D $\leftrightarrow$ K) interchanges.
Because the D, E, R, and K residues engage in nonelectrostatic as well as
electrostatic interactions, sequence-specific nonelectrostatic interactions do
contribute to $\Delta(T_{\rm cr})_{\rm WT-CS}$ and this effect clearly depends
on the interaction scheme. Now, the energetic effects of the (E
$\leftrightarrow$ R), (D $\leftrightarrow$ R), (E $\leftrightarrow$ K), and (D
$\leftrightarrow$ K) swaps may be characterized by ${\cal E}\Delta\lambda_{\rm
E-R}/2$, ${\cal E}\Delta\lambda_{\rm D-R}/2$, ${\cal E}\Delta\lambda_{\rm
E-K}/2$, and ${\cal E}\Delta\lambda_{\rm D-K}/2$; but these quantities can
have different magnitudes and even different signs for different interaction
schemes, as exemplified by the ${\cal E}\Delta\lambda_{\rm E-R}/2$, ${\cal
E}\Delta\lambda_{\rm D-R}/2$, ${\cal E}\Delta\lambda_{\rm E-K}/2$, and ${\cal
E}\Delta\lambda_{\rm D-K}/2$ values of $0.0459$, $0.0378$, $-0.0055$, and
$-0.0136$ kcal mol-1, respectively, for HPS, the corresponding values of
$-0.0558824$, $-0.0264705$, $-0.0382354$, and $-0.0088235$ kcal mol-1 for
Urry, and $0.018596$, $0.0065$, $-0.004485$, and $-0.016581$ for FB.
Ramifications of these differences on the LLPS properties of model CS Ddx4 IDR
remain to be further explored.

Coarse-Grained Explicit-Chain MD for Hydrophobic-Polar Sequences. Fig. 12
provides phase behaviors computed by explicit-chain MD simulation for the
three LS sequences studied above by FTS. The rank ordering of LLPS propensity
LS3 $>$ LS2 $\gtrsim$ LS1 observed in Fig. 12a is consistent with that deduced
from FTS correlation functions and PMFs in Fig. 8, attesting to the
effectiveness of the FTS formulation developed here for sequence-specific
effects of spatially short-range interactions on LLPS of IDPs. There are minor
differences between the FTS and MD results, in that the FTS correlation
functions and PMFs suggest that the LLPS propensities of LS1 and LS2 are
practically identical considering numerical uncertainties, the MD phase
diagrams show that LS2 has a slightly higher LLPS propensity than LS1. Such
minor differences are not unexpected, however, given the differences discussed
above between the FTS and MD models. In both our FTS and MD models, the more
blocky LS3 sequence clearly exhibits its superior LLPS propensity.

Fig. 12: Sequence-dependent phase behaviors of hydrophobic-polar sequences in
an explicit-chain model. (a) Phase diagrams (coexistence curves) of the three
leucine-serine sequences LS1, LS2, and LS3 in Fig. 8a are computed by coarse-
grained MD using the KH interaction schemedignon18 and the protocol described
in Models and Methods. Simulation data points and extrapolated critical
pointsdignon18 ; SumanPNAS ; panag2017 are shown as filled circles, connecting
lines are a guide for the eye. (b--d) Snapshots taken during the simulation
showing a part of the simulation box for 1,000 copies of (b) LS1 at model
temperature $T=1000$ K, (c) LS2 at model $T=1000$ K, and (d) LS3 at model
$T=1200$ K. Beads representing leucine and serine are shown, respectively, in
orange and cyan in these snapshots. Significant number of polymers are seen in
the dilute phase of each of these systems.

The critical temperatures of the LS sequences shown in Fig. 12a are $3$\--$4$
times higher than those of the Ddx4 IDRs in Fig. 11, exhibiting stable
droplets at model temperatures as high as $1000$ K (Fig. 12b--d). This
observation is notable because LLPS propensities tend to increase with chain
lengthlinPRL but the LS sequences have a much shorter chain length ($N=20$)
than that of the Ddx4 IDRs ($N=236$). A likely reason for the high LLPS
propensities of the LS sequences is the strongly favorable leucine-leucine
interaction (strongest in the KH scheme) and favorable leucine-serine
interactions in the model. The snapshot in Fig. 12d for LS3 shows a certain
degree of separate clustering of the leucine residues and of the serine
residues, but not to the extent displayed by the micellar organization of the
corresponding hydrophobic-polar sequence in ref. (126), probably because the
leucine-serine interaction is favorable in our KH interaction scheme (and thus
conducive to leucine-serine mixing) but the interaction between a hydrophobic
and a polar bead is repulsive in ref. (126).

A sequence hydropathy decoration (SHD) parameter, defined as SHD
$=N^{-1}\sum_{\alpha<\beta}^{N}[\lambda(r_{\alpha})+\lambda(r_{\beta})]|\alpha-\beta|^{-1}$,
was proposed recently as a predictor of properties of single-chain
conformational ensembles of heteropolymers.zhengHP For applications to
interaction schemes such as KH and Mpipi that have 210 contact energies
instead of a 20-value hydrophobicity/hydropathy scale, we consider a natural
generalization of the above formula:

${\rm SHD}\rightarrow{\rm SHD}^{\rm
hh}\equiv-N^{-1}\sum^{N}_{\alpha<\beta}\varepsilon_{r_{\alpha},r_{\beta}}|\alpha-\beta|^{-1}\;,$
(62)

where $\varepsilon_{r,r^{\prime}}$ for KH is given by the entries in Table S3
of ref. (40). For the LS1, LS2, and LS3 sequences we consider here in the KH
interaction scheme, ${\rm SHD}^{\rm hh}=1.12$, $1.05$, and $1.17$,
respectively. As the ${\rm SHD}^{\rm hh}$ value of LS3 with higher LLPS
propensity is larger than those of LS1 and LS2, this result suggests that SHD
or ${\rm SHD}^{\rm hh}$ may be used as a predictor for LLPS as well. Further
effort will be needed to examine whether the correlation between SHD or ${\rm
SHD}^{\rm hh}$ and LLPS propensity of heteropolymers with spatially short-
range interactions is as strong as that between sequence charge decoration
(SCD)kings2015 and LLPS propensity of polyampholytes.lin2017 ; suman1 ; suman2

CONCLUSIONS

In summary, we have developed a field-theoretic formulation for modeling
sequence-specific biomolecular phase separation. Our theory offers a coarse-
grained account of short-spatial-range $\pi$-related and hydrophobic
interactions as well as long-spatial-range Coulomb interactions. In
conjunction with RPA for electrostatics, a mean-field approximation for
spatially short-range interactions derived from the general theory is useful
as a computationally efficient tool for studying phase separation of
intrinsically disordered proteins, as exemplified by the application to the
Ddx4 IDRs described here. Full sequence effects of spatially short-range
interactions can be studied using FTS, illustrated here by the sequence-
dependent phase properties of three different hydrophobic-polar sequences with
the same hydrophobic/polar composition. Initial success in applying these
approaches to several different interaction schemes for modeling biomolecular
LLPS are confirmed by coarse-grained explicit-chain molecular dynamics
simulations. We have included only Yukawa potentials in our theory for this
initial effort. Because ostensibly small changes in coarse-grained interaction
potentials can lead to fundamental variations in conformational
properties,cheung2002 ; liu2005 ; kaya2013 future effort should aim to extend
our formulation to other functional forms for spatially short-range
potentials.ottinger2021 Techniques should also be further developed to study
longer sequences of biological IDRs by FTS. Much exciting work lies ahead
under the present theoretical framework.

Acknowledgements. We thank Yi-Hsuan Lin for helpful discussions. Financial
support for this work was provided by Canadian Institutes of Health Research
grant NJT-155930 and Natural Sciences and Engineering Research Council of
Canada Discovery grant RGPIN-2018-04351 to H.S.C. We are grateful for the
computational resources provided generously by Compute/Calcul Canada.

The authors declare no conflict of interest.

References

## References

  * (1) Banani, S. F.; Lee, H. O.; Hyman, A. A.; Rosen, M. K. Biomolecular condensates: organizers of cellular biochemistry. Nat. Rev. Mol. Cell Biol. 2017, 18, 285--298.
  * (2) Shin, Y.; Brangwynne, C.P. Liquid phase condensation in cell physiology and disease. Science 2017, 357, eaaf4382.
  * (3) Lyon, A. S.; Peeples, W. B.; Rosen, M. K. A framework for understanding the functions of biomolecular condensates across scales. Nat. Rev. Mol. Cell Biol. 2021, 22, 215--235.
  * (4) Jain, S.; Wheeler, J. R.; Walters, R. W.; Agrawal, A.; Barsic, A.; Parker, R. ATPase-modulated stress granules contain a diverse proteome and substructure. Cell 2016, 164, 487--498.
  * (5) Wilson, E. B. The structure of protoplasm. Science 1899, 10, 33--45.
  * (6) Brangwynne, C. P.; Eckmann, C. R.; Courson, D. S.; Rybarska, A.; Hoege, C.; Gharakhani, J.; Jülicher, F.; Hyman, A. A. Germline P granules are liquid droplets that localize by controlled dissolution/condensation. Science 2009, 324, 1729--1732.
  * (7) Li, P.; Banjade, S.; Cheng, H. C.; Kim, S.; Chen, B.; Guo, L.; Llaguno, M.; Hollingsworth, J. V.; King, D. S.; Banani, S. F.; Russo, P. S.; Jiang, Q.-X.; Nixon, B. T.; Rosen, M. K. Phase transitions in the assembly of multivalent signalling proteins. Nature 2012, 483, 336--340.
  * (8) Kato, M.; Han, T. W.; Xie, S.; Shi, K.; Du, X.; Wu, L. C.; Mirzaei, H.; Goldsmith, E. J.; Longgood, J.; Pei, J.; Grishin, N. V.; Franz, D. E.; Schneider, J. W.; Chen, S.; Li, L.; Sawaya, M. R.; Eisenberg, D.; Tycko, R.; McKnight, S. L. Cell-free formation of RNA granules: low complexity sequence domains form dynamic fibers within hydrogels. Cell 2012, 149, 753--767.
  * (9) Nott, T. J.; Petsalaki, E.; Farber, P.; Jervis, D.; Fussner, E.; Plochowietz, A.; Craggs, T. D.; Bazett-Jones, D. P.; Pawson, T.; Forman-Kay, J. D.; Baldwin, A. J. Phase transition of a disordered nuage protein generates environmentally responsive membraneless organelles. Mol. Cell 2015, 57, 936--947.
  * (10) Molliex, A.; Temirov, J.; Lee, J.; Coughlin, M.; Kanagaraj, A. P.; Kim, H. J.; Mittag, T.; Taylor, J. P. Phase separation by low complexity domains promotes stress granule assembly and drives pathological fibrillization. Cell 2015, 163, 123--133.
  * (11) Elbaum-Garfinkle, S.; Kim, Y.; Szczepaniak, K.; Chen, C. C.-H.; Eckmann, C. R.; Myong, S.; Brangwynne, C. P. The disordered P granule protein LAF-1 drives phase separation into droplets with tunable viscosity and dynamics. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, 7189--7194.
  * (12) Lin, Y.; Protter, D. S. W.; Rosen, M. K.; Parker, R. Formation and maturation of phase-separated liquid droplets by RNA-binding proteins. Mol. Cell 2015, 60, 208--219.
  * (13) Bergeron-Sandoval, L.P.; Safaee, N.; Michnick, S. W. Mechanisms and consequences of macromolecular phase separation. Cell 2016, 165, 1067--1079.
  * (14) Li, X.-H.; Chavali, P. L.; Pancsa, R.; Chavali, S.; Babu, M. M. Function and regulation of phase-separated biological condensates. Biochemistry 2018, 57, 2452--2461.
  * (15) Boeynaems, S.; Alberti, S.; Fawzi, N. L.; Mittag, T.; Polymenidou, M.; Rousseau, F.; Schymkowitz, J.; Shorter, J.; Wolozin, B.; Van Den Bosch, L.; Tompa, P.; Fuxreiter, M. Protein phase separation: a new phase in cell biology. Trends Cell Biol. 2018, 28, 420--435.
  * (16) Cinar, H.; Fetahaj, Z.; Cinar, S.; Vernon, R. M.; Chan, H. S.; Winter, R. Temperature, hydrostatic pressure, and osmolyte effects on liquid-liquid phase separation in protein condensates: Physical chemistry and biological implications. Chem. Eur. J. 2019, 57, 13049--13069.
  * (17) Gomes, E., and Shorter, J. The molecular language of membraneless organelles. J. Biol. Chem. 2019, 294, 7115--7127.
  * (18) Burke, K. A.; Janke, A. M.; Fhine, C. L.; Fawzi, N. L. Residue-by-residue view of in vitro FUS granules that bind the C-terminal domain of RNA polymerase II. Mol. Cell 2015, 60, 231--241.
  * (19) Brady, J. P.; Farber, P. J.; Sekhar, A.; Lin, Y.-H.; Huang, R.; Bah, A.; Nott, T. J.; Chan, H. S.; Baldwin, A. J.; Forman-Kay, J. D.; Kay, L. E. Structural and hydrodynamic properties of an intrinsically disordered region of a germ cell-specific protein on phase separation. Proc. Natl. Acad. Sci. U.S.A. 2017, 114, E8194--E8203.
  * (20) Kato, M.; McKnight, S. L. A solid-state conceptualization of information transfer from gene to message to protein. Annu. Rev. Biochem., 2018, 87, 351-390.
  * (21) Feric, M.; Vaidya, N.; Harmon, T. S.; Mitrea, D. M.; Zhu, L.; Richardson, T. M.; Kriwacki, R. W.; Pappu, R. V.; Brangwynne, C. P. Coexisting liquid phases underlie nucleolar subcompartments. Cell 2016, 165, 1686--1697.
  * (22) Lin, Y. -H.; Brady, J. P.; Forman-Kay J. D.; Chan, H. S. Charge pattern matching as a ‘fuzzy’ mode of molecular recognition for the functional phase separations of intrinsically disordered proteins. New J. Phys. 2017, 19, 115003.
  * (23) Jülicher, F.; Grill, S. W.; Salbreux, G. Hydrodynamic theory of active matter. Rep. Prog. Phys. 2018, 81, 076601.
  * (24) Bertrand, T.; Lee, C. F. Diversity of phase transitions and phase separations in active fluids. Phys. Rev. Res. 2022, 4, L022046.
  * (25) Rosowski, K. A.; Sai, T.; Vidal-Henriquez, E.; Zwicker, D.; Style, R. W.; Dufresne, E. R. Elastic ripening and inhibition of liquid–liquid phase separation. Nat. Phys. 2020, 16, 422--425.
  * (26) Zhou, H.-X. Viscoelasticity of biomolecular condensates conforms to the Jeffreys model. J. Chem. Phys. 2021, 154, 041103.
  * (27) Sanchez-Burgos, I.; Joseph, J. A.; Collepardo-Guevara, R.; Espinosa, J. R. Size conservation emerges spontaneously in biomolecular condensates formed by scaffolds and surfactant clients. Sci. Rep. 2021, 11, 15241.
  * (28) Harmon, T. S.; Holehouse, A. S.; Rosen, M. K.; Pappu, R. V. Intrinsically disordered linkers determine the interplay between phase separation and gelation in multivalent proteins. eLife 2017, 6, e30294.
  * (29) Lin, Y.-H.; Forman-Kay, J. D.; Chan, H. S. Theories for sequence-dependent phase behaviors of biomolecular condensates. Biochemistry 2018, 57, 2499--2508.
  * (30) McSwiggen, D. T.; Mir, M.; Darzacq, X.; Tjian, R. Evaluating phase separation in live cells: Diagnosis, caveats, and functional consequences. Genes Dev. 2019, 33, 1619--1634.
  * (31) Musacchio, A. On the role of phase separation in the biogenesis of membraneless compartments. EMBO J. 2022, 41, e109952.
  * (32) Mittag, T.; Pappu, R. V. A conceptual framework for understanding phase separation and addressing open questions and challenges. Mol. Cell 2022, 82, 2201--2214.
  * (33) Lin, Y.-H.; Wu, H.; Jia, B.; Zhang, M.; Chan, H. S. Assembly of model postsynaptic densities involves interactions auxiliary to stoichiometric binding. Biophys. J. 2022, 121, 157--171.
  * (34) Hyman, A. A.; Weber, C. A.; Jülicher, F. Liquid-liquid phase separation in biology. Annu. Rev. Cell Dev. Biol. 2014, 30, 39--58.
  * (35) Brangwynne, C. P.; Tompa, P.; Pappu, R. V. Polymer physics of intracellular phase transitions. Nat. Phys. 2015, 11, 899--904.
  * (36) Lin, Y.-H.; Forman-Kay, J. D.; Chan, H. S. Sequence-specific polyampholyte phase separation in membraneless organelles. Phys. Rev. Lett., 2016, 117, 178101.
  * (37) Vovk, A.; Gu C.; Opferman, M. G.; Kapinos, L. E.; Lim, R. Y. H.; Coalson, R. D.; Jasnow, D.; Zilman, A. Simple biophysics underpins collective conformations of the intrinsically disordered proteins of the Nuclear Pore Complex. eLife 2016, 5, e10785.
  * (38) Lin, Y.-H.; Song, J.; Forman-Kay, J. D.; Chan, H. S. Random-phase-approximation theory for sequence-dependent, biologically functional liquid-liquid phase separation of intrinsically disordered proteins. J. Mol. Liq., 2017, 228, 176--193.
  * (39) Lin, Y.-H.; Chan, H. S. Phase separation and single-chain compactness of charged disordered proteins are strongly correlated. Biophys, J. 2017, 112, 2043--2046.
  * (40) Dignon, G. L.; Zheng, W.; Kim, Y. C.; Best, R. B.; Mittal, J. Sequence determinants of protein phase behavior from a coarse-grained model. PLoS Comput. Biol. 2018, 14, e1005941.
  * (41) Das, S.; Eisen, A.; Lin, Y.-H.; Chan, H. S. A lattice model of charge-pattern-dependent polyampholyte phase separation. J. Phys. Chem. B 2018, 122, 5418--5431.
  * (42) Dignon, G. L.; Zheng, W.; Best, R. B.; Kim, Y. C.; Mittal, J. Relation between single-molecule properties and phase behavior of intrinsically disordered proteins. Proc. Natl. Acad. Sci. U.S.A. 2018, 115, 9929--9934.
  * (43) Das, S.; Amin, A. N.; Lin, Y.-H.; Chan, H. S. Coarse-grained residue-based models of disordered protein condensates: Utility and limitations of simple charge pattern parameters. Phys. Chem. Chem. Phys. 2018, 20, 28558--28574.
  * (44) Robichaud,N. A. S.; Saika-Voivod, I.; Wallin, S. Phase behavior of blocky charge lattice polymers: Crystals, liquids, sheets, filaments, and clusters. Phys. Rev. E. 2019, 100, 052404.
  * (45) Lin, Y.; McCarty, J.; Rauch, J. N.; Delaney, K. T.; Kosik, K. S.; Fredrickson, G. H.; Shea, J.-E.; Han, S. Narrow equilibrium window for complex coacervation of tau and RNA under cellular conditions eLife 2019, 8, e42571.
  * (46) McCarty, J.; Delaney, K. T.; Danielsen, S. P. O.; Fredrickson, G. H.; Shea, J.-E. Complete phase diagram for liquid-liquid phase separation of intrinsically disordered proteins. J. Phys. Chem. Lett. 2019, 10, 1644--1652.
  * (47) Leopold, P. E.; Montal, M.; Onuchic, J. N. Protein folding funnels: A kinetic approach to the sequence-structure relationship. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 8721--8725.
  * (48) Wolynes, P. G.; Onuchic, J. N.; Thirumalai, D. Navigating the folding routes. Science 1995, 267, 1619--1620.
  * (49) Dill, K. A.; Chan, H. S. From Levinthal to pathways to funnels. Nat. Struct. Biol. 1997, 4, 10--19.
  * (50) Klosin, A.; Oltsch, F.; Harmon, T.; Honigmann, A.; Jülicher, F.; Hyman, A. A.; Zechner, C. Phase separation provides a mechanism to reduce noise in cells. Science 2020, 367, 464--468.
  * (51) Deviri, D.; Safran, S. A. Physical theory of biological noise buffering by multicomponent phase separation. Proc. Natl. Acad. Sci. U.S.A. 2021, 118, e2100099118.
  * (52) Cinar, H.; Oliva, R.; Lin, Y.-H.; Chen, X.; Zhang, M.; Chan, H. S.; Winter, R. Pressure sensitivity of SynGAP/PSD-95 condensates as a model for postsynaptic densities and its biophysical and neurological ramifications. Chem. Eur. J. 2020, 26, 11024--11031.
  * (53) Cinar, H.; Oliva, R.; Wu, H.; Zhang, M.; Chan, H. S.; Winter, R. Effects of cosolvents and crowding agents on the stability and phase transition kinetics of the SynGAP/PSD-95 condensate model of postsynaptic densities. J. Phys. Chem. B 2022, 126, 1734--1741.
  * (54) Tsang, B.; Pritišanac, I.; Scherer, S. W.; Moses, A. M.; Forman-Kay, J. D. Phase separation as a missing mechanism for interpretation of disease mutations. Cell 2020, 183, 1742--1756.
  * (55) Flory, P. J. (1953) Principles of Polymer Chemistry; Cornell University Press:Ithaca, New York, 1953.
  * (56) Overbeek, J. T. G.; Voorn, M. J. Phase separation in polyelectrolyte solutions. Theory of complex coacervation. J. Cell. Comp. Physiol. 1957, 49, 7--26.
  * (57) Ermoshkin, A. V.; Olvera de la Cruz, M. A modified random phase approximation of polyelectrolyte solutions. Macromolecules 2003, 36, 7824--7832.
  * (58) Fredrickson, G. H. The Equilibrium Theory Of Inhomogeneous Polymers; Oxford University Press Inc.:New York, 2006.
  * (59) Danielsen, S. P. O.; McCarty, J.; Shea, J.-E.; Delaney, K. T.; Fredrickson, G. H. Molecular design of self-coacervation phenomena in block polyampholytes. Proc. Natl. Acad. Sci. U.S.A. 2019, 116, 8224--8232.
  * (60) Danielsen, S. P. O.; McCarty, J.; Shea, J.-E.; Delaney, K. T.; Fredrickson, G. H. Small ion effects on self-coacervation phenomena in block polyampholytes. J. Chem. Phys. 2019, 151, 034904.
  * (61) Nilsson, D.; Irbäck, A. Finite-size scaling analysis of protein droplet formation. Phys. Rev. E, 2020, 101, 022413.
  * (62) Pal, T.; Wessén, J.; Das, S.; Chan, H. S. Subcompartmentalization of polyampholyte species in organelle-like condensates is promoted by charge-pattern mismatch and strong excluded-volume interaction. Phys. Rev. E 2021, 103, 042406.
  * (63) Nilsson, D.; Irbäck, A. Finite-size shifts in simulated protein droplet phase diagrams. J. Chem. Phys. 2021, 154, 235101.
  * (64) Nilsson, D.; Bozorg, B.; Mohanty, S.; Söderberg, B.; Irbäck, A. Limitations of field-theory simulation for exploring phase separation: The role of repulsion in a lattice protein model. J. Chem. Phys. 2022, 156, 015101.
  * (65) Choi, J.-M.; Dar, F.; Pappu, R. V. LASSI: A lattice model for simulating phase transitions of multivalent proteins. PLoS Comput. Biol. 2019, 15, e1007028.
  * (66) Das, S.; Lin, Y.-H.; Vernon, R. M.; Forman-Kay, J. D.; Chan, H. S. Comparative roles of charge, $\pi$, and hydrophobic interactions in sequence-dependent phase separation of intrinsically disordered proteins. Proc. Natl. Acad. Sci. U.S.A. 2020, 117, 28795--28805.