Analytical Formulation and Field-Theoretic Simulation of Sequence-Specific Phase Separation of Proteinlike Heteropolymers with Short- and Long-Spatial-Range Interactions

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 sequence-pattern 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 a set of different model hydrophobic-polar sequences sharing the same composition. Effectiveness of the methodology is 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. [1][2][3] 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 ∼ 300 protein species identified from proteomic analyses of yeast and mammalian stress granules using mass spectrometry. 4 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, 5 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. 6 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. [7][8][9][10][11][12][13][14][15][16][17] One novel observation from these investigations is that multivalent, relatively weak interactions of IDPs and IDRs-typically weaker than binding of folded domains-often play a major role in the formation of biomolecular condensates. [8][9][10]12 These interactions among largely disordered conformations of IDPs and IDRs are transient and presumably nonstoichiometric, 18,19 though in some instances they entail sampling disordered chain configurations with labile but nonetheless specific fibril-like local structures 8,20 . Despite the stochastic nature of these IDP/IDR interactions, they are capable of contributing to a "fuzzy", sequence-dependent mechanism of molecular recognition. 21,22 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; 23,24 and formulations beyond simple LLPS theory are needed to describe the size, viscoelastic, and other materials properties of biomolecular condensates. [25][26][27] In addition to physicochemical processes akin to LLPS, gelation and percolation 28,29 as well as structure-specific stoichiometric interactions 30,31 are expected to also contribute significantly to the formation of biomolecular condensates. It is likely that these processes are coupled 32 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. 33 Nonetheless, despite these complexities, since LLPS is integral to-though not exclusively responsible for-biomolecular condensate formation, experimental [6][7][8][9][10][11][12][13][14][15][16][17] and theoretical/computational 22,[34][35][36][38][39][40][41][42][43][44][45][46] 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 scienceincluding, as a befitting example for this Special Issue, a seminal lattice model study by Onuchic and coworkers 47 that inspired the funnel landscape picture of protein folding, 48,49 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, 22,29,50,51 novel clues to neurological effects of hydrostatic pressure, 52,53 and a likely link between neurological disease-causing mutations and bioinformatics-inferred LLPS propensities, 54 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. 29 These modeling approaches include-but certainly not limited to-mean-field Flory-Huggins (FH) 55 and Overbeek-Voorn (OV) 56 theories of polymer solutions, random phase approximation (RPA) theories to address effects of sequence charge patterns on LLPS of polyampholytes as model IDPs/IDRs, 36,38,57 field-theoretic simulation (FTS) 58 to provide an improved account of charge and matter density fluctuations beyond RPA for polyampholyte LLPS 45,46,59,60 -though FTS is still limited by finite-size effects and in its treatment of excluded volume 61-64 -as well as explicit-chain simulations of lattice, 21,41,44,65 continuum coarse-grained, 40,42,43,[66][67][68] and atomistic 69,70 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; [71][72][73] and LLPS of single or multiple biomolecular species, especially those involving folded proteins and folded domains, have also been modeled by patchy particles [74][75][76] using simulations 75,76 as well as analytical formalisms based on Wertheim's thermodynamic perturbation theory. 74,77 As expected, the required computation increases with structural and energetic details that a model seeks to capture. In this context, it is noteworthy that even basic FH theory, which requires minimal numerical effort 78 yet has recently been made even more tractable by a novel self-consistent solution, 79 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, 22,29,51 to ascertain the extent of void-volume contributions to the hydrostatic pressure dependence of LLPS, 52 and to rationalize experimental data on the impact of aromatic valence on IDP LLPS. 80 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, [40][41][42][43][44]66,69,70,[81][82][83][84][85] 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, 86 and RPA augmented by Kuhn-length renormalization (rG-RPA) 87 for a more physical account of charge density and its fluctuations, 88 have focused only on spatially long-range Coulomb interactions entailed by the heteropolymeric sequence patterns of electric charges 86,89,90 and their effects on LLPS, 36,38,39,45,46,59,60,62,78,87,91,92 single-chain polyampholyte/IDP conformational properties, 39,88,90,[93][94][95][96] and the interaction of a pair of polyampholyte chains. 81 These theoretical treatments of spatially long-range electrostatic interactions have led to many important physical insights; but these theories by themselves do not address the full sequence-dependent effects of spatially short-range, contact-like hydrophobic 69 and π-related 97,98 interactions, which are integral parts of the physical driving forces for biomolecular LLPS 40,66,69,[98][99][100][101] as well as IDP conformational properties 42,102,103 and interactions. 97 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. 36,38 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 free energy, signalling the instability of the homogeneous phase about which the RPA expansion is performed. 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, properly handle the 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 energies 104-107 from either the values for the 20 amino acids on a hydrophobicity (hydropathy) scale 40,[108][109][110][111][112][113] or the 20 × 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. 40,[114][115][116][117][118] This approach is computationally tractable because matrices for physical pairwise interactions among amino acid residues tend to possess a dominant eigenvector. [104][105][106]119,120 Thus, only a small number of types of nonelectric charge ( 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, 121-123 the size-and length-scale dependence of hydrophobic interactions [124][125][126] and their deviations from pairwise additivity 127-129 (as manifested, e.g., by the fact that pairwise interactions alone are insufficient to account for protein folding cooperativity 127,130,131 ), we apply our new formulation to several different hydrophobicity scales 40,111,112 and pairwise interaction schemes 40,117 to compare their predicted LLPS behaviors for the wildtype and variant IDRs of the DEAD-box RNA helicase Ddx4 9,19 and two-letter (hydrophobic-polar) model sequences. 132,133 The theoretical predictions are compared against experimental trend and corresponding predictions from coarse-grained explicit-chain molecular dynamics 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 × 20 matrix ε r,r of contact energies for the pairwise interactions between residue types r and r , this matrix can be written in the diagonalized form because ε r,r is symmetric (ε r,r = ε r ,r ), with {λ a } 20 a=1 being the eigenvalues of ε r,r , 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., |λ a | ≥ |λ a+1 |.
The spectral form in Eq. (1) implies that we can decompose the full ε r,r into a set of interactions, with overall strengths λ a and relative residue-dependent strengths given by product of nonelectric charges q (a) r q (a) r . 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. 104,105,114,115 Although all twenty eigenvalues are needed in Eq. (1) to recover the full ε r,r exactly, one can often achieve a very good numerical agreement by including only the first few most dominant λ a , i.e. by truncating the sum at someã < 20, because of the pattern of hydrophobic-like interactions among amino acid residues. [104][105][106]119,120 For our present purpose, we consider two interaction matrices proposed recently for the study of IDP LLPS, namely, ε KH r,r obtained through a shifted MJ matrix, 40,115 and ε Mpipi r,r , 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 IDPs 117 (Fig. 1). The accuracies of truncated eigenvalue decompositions of ε KH r,r and ε Mpipi r,r [Eq.
(2), withã = 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 shows significant deviations from ε KH r,r or ε Mpipi r,r when only one eigenvalue is used (ã = 1). At the same time, it indicates that eigenvalue decomposition is reasonably accurate forã ≥ 2. Interestingly, forã = 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, 104,105,120 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 ≥ 2. Intuitively, for KH, the nonelectric charge for a = 2 corresponds to a measure of polarity, with most polar residues   1) and (2) where the summations over eigenvalues are truncated atã. 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 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 tentaliz-ing, 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 representative 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) ε HPS r,r , which is derived in ref.
whereλ r are twenty residue-specific hydrophobicity/hydropathy values [not to be confused with the eigenvalues λ a in Eq. (1)] and ∆ is an overall shift. Matrices in the form of Eq. (3) have at most two nonzero eigenvalues, given by λ 1,2 = 1 2 rλ r + ∆ ± 20 r (λ r + ∆) 2 . 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.
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 × 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 ε KH r,r and ε Mpipi r,r is assessed by considering the quantity which is the difference between the truncated summation and the original ε r,r . Results for KH and Mpipi atã = 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ã = 3 with the only exception of the cysteine-cysteine (C-C) entry, which is noticeably overestimated at allã ≤ 3. This exception is in line with the fact that the covalent C-C disulfide bond   is physically different from all the 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ã = 3 is reasonably accurate though it overestimates interactions involving the electrically negative aspartic and glutamic acid residues.
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 ε r,r . Relative to the field theory we consider previously for polyampholytes, 62,78,91,92 the present extended field theory contains a additional fields, i.e., one additional field ϕ a (r) for every eigenvalue λ a included, wherein the field ϕ a (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 ε r,r 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. Statistical Field Theory. Here we consider a system of volume V with n p protein chains, each consisting of N consecutive amino acid residues of types r α (with α = 1, . . . , N ), and denote the position of residue α on chain i (= 1, . . . , n p ) as R i,α . The system contains also n + unit-positive and n − unit-negative ions with positions r ±,i , and n w neutral solvent particles with positions r w,i . We let ρ p,±,w ≡ n p,±,w /V denote the bulk number densities of the polymers, ions and solvent in the system, and ρ b = N n p /V = N ρ p the bulk number density of polymer beads. The canonical partition function of this system is given by

Mpipi
where are position integration measures for polymer beads, ions and solvents, respectively. The microscopic Hamiltonian in units of k B T , where k B is Boltzmann's constant and T is absolute temperature, written asĤ contains terms accounting for chain connectivity (Ĥ 0 ), soft compressibility (Ĥ c ), electrostatic interactions (Ĥ e ) and spatially short-range hydrophobic-like interactions-which may include van der Waals, π-related, and other forms of short-spatial-range interactions (Ĥ h ). The chain connectivity termĤ 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 penalizes deviations of the total densitŷ from a reference density ρ 0 , and thus serves to capture an essential aspect of excluded volume effects. In the above equation,ρ i (r) with i = b, +, −, w (referring to polymer bead, cation, anion, and solvent particle, respectively) is the number density of species i at position r and the v i s are factors that we use to model the relative volume of the different particles. The number densities arê where individual particles are modeled as Gaussian distributions Γ(r) = e −r 2 /2a 2 s /(2πa 2 s ) 3/2 , which serve to regularize infinities arising from particle contact interactions and self interactions. 135,136 Here the smearing parameter a s is equal to 1/ √ 3 of the standard deviation of r in the Γ(r) distribution (i.e., a s = |r| 2 /3). We refer to γ in Eq. (9)-which has units of density-as the compressibility, and note that in the limit of γ → 0, we obtain the standard incompressibility condition e −Ĥc ∝ δ[ρ(r) − ρ 0 ] where the functional δ-function enforcesρ(r) = ρ 0 at all spatial positions r.
Electrostatic interactions are included inĤ e through the standard pairwise Coulomb potential,Ĥ whereĉ is the density of electric charge at r and σ i,α is the electric charge of the αth residue on the ith chain in units of the proton charge e. Because the present study is limited to systems of chains with identical sequence, i.e., σ i,α = σ α is independent of i, the subscript i will be dropped from the symbol σ for electric charge hereafter. The strength of electrostatic interactions is controlled by the Bjerrum length l B = e 2 /4π 0 r k B T that involves vacuum permittivity 0 and a constant background relative permittivity r . In this work, we neglect any concentration dependence of r as this has been shown to amount to only small to moderate effects on LLPS propensity. 91 Nonetheless, if desired, it would be straightforward to implement a concentration dependent permittivity in our model through the approach described in ref. 91. The final term in the Hamiltonian in Eq. (7),Ĥ h , contains the pairwise residue-specific spatially short-range interactions, which we express aŝ where accounts for the residue-residue interaction in the presence of Gaussian smearing (note that . Leaving the functional form of the interaction potential V h (|r|) (in units of k B T ) unspecified for the moment, we have in the above expression assumed that all residue dependence is accounted for by the overall multiplicative factors ε rα,r β from a given interaction matrix. The spectral form of ε r,r in Eq. (1) can be used to writeĤ h in Eq. (14) asĤ whereĥ is the density of nonelectric charge for the ath eigenvalue λ a at position r.
The Hamiltonian termsĤ c,e,h for non-bonded interactions are quadratic in densities and can thus be decoupled using standard Hubbard-Stratonovich transformations, 58 turning the partition function in Eq. (5) into that of a statistical field theory. Special care needs to be given toĤ h wherein terms with λ a > 0 and λ a < 0 need to be linearized with fields integrated on the imaginary and real axis, respectively. To keep track of the sign of eigenvalue λ a , we introduce the variable ξ a , defined as where 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 where the field Hamiltonian (in units of k B T ) In the field theory representation of the system, η(r) and ψ(r) are conjugate fields to the densityρ(r) and the electric charge densityĉ(r), respectively, while ϕ a (r) is conjugate to the nonelectric charge densityĥ a (r) associated with eigenvalue λ a in the spectral decomposition of the interaction matrix. The single-molecule partition functions Q p,±,w [η,ψ, {φ a }] depend on smeared fieldsφ(r) = Γ φ(r) ≡ dr Γ(r − r )φ(r ) (for φ = η, ψ, ϕ a ) via the relations for ions and solvents, respectively, and for polymer chains, where N ≡ V (2πb 2 /3) 3(N −1)/2 normalizes Q p to unity at zero field values, and In Eq. (20), V h −1 represents a differential operator satisfying and thus V h −1 depends on the functional form of V h (|r|). Since our primary goal here is formulation development, for simplicity we consider only temperature-independent forms of the spatially short-range interaction potential V h (|r|). Physically, however, solvent-mediated interactions of amino acid residues can be temperature dependent when viewed in an implicit-solvent perspective, 29 especially for the hydrophobic effects among nonpolar residues, 49 leading, e.g., to cold denaturation 137 and LLPS driven by increasing temperature rather than decreasing temperature. 16 When necessary, these effects can be readily incorporated into our formulation by introducing temperature-dependent potentials, 49,137,138 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 η and ϕ a , while electrostatic effects first show up at the one-loop level in ψ which are accounted for by RPA.
The spatially homogeneous MFT solution is obtained by first letting the fields η, ψ and ϕ a in the expression for the field Hamiltonian H in Eq. (20) be position-independent (hence smeared fieldsφ = φ for φ = η, ψ, ϕ a are also position-independent) and then requiring that the first derivatives of the resulting H with respect to the fields are zero. With positionindependent fields, the single-molecule partition functions Q i , defined in Eqs.
where σ p = N α=1 σ α is the net charge of a polymer and q (a) = N α=1 q (a) rα is the total nonelectric charge associated with the ath eigenvalue of one polymer chain. The above relations for Q i lead to the following expression for H with homogeneous fields: In Eq. (27), of the spatially short-range potential V h (|r|). This factor can be derived by considering the Then, by noting that the Fourier transform of the position-independent ϕ a is proportional to δ(k), it is straightforward to show that the dr integration with position-independent ϕ a produces theV h (0) term.
The coefficient of iψ in Eq. (27) is the total electric charge in the system which must be set to zero, n p σ p + n + − n − = 0. The charge-conjugate field ψ thus vanishes in Eq. (27) with no contribution to the MFT solution. For the other fields, solving ∂H/∂η = ∂H/∂ϕ a = 0 forη andφ a leads toη where ∆ρ ≡ ν b N ρ p +ν + ρ + +ν − ρ − +ν w ρ w −ρ 0 is the deviation of the total bulk density from ρ 0 , and q (a) ≡ N α=1 q (a) rα is the total nonelectric charge associated with the ath eigenvalue λ a of one polymer chain. The same expressions forη andφ a [Eq. (28)] may alternatively be obtained by requiring vanishing first functional derivatives of H with respect to η(r) and ϕ a (r) and then seeking position-independent solutions η(r) →η, ψ(r) →ψ and ϕ a (r) →φ a of the resulting δH/δη(r) = δH/δψ(r) = δH/δϕ a (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 volumeh ≡ H/V , given bȳ In the approximate analytical approach taken in this work, we set η(r) and ϕ a (r) to their MFT valuesη andφ a , while a spatially varying ψ(r) is kept to O(ψ 2 ) in H. The field Hamiltonian in this approximation is therefore The factor g cc (k), given by 91 comes from the quadratic expansion of the single polymer partition function Q p . In our approximation, the functional integrals in the partition function Z amounts to where Vh is the field Hamiltonian evaluated at homogeneous field values in accordance with Eq. (29).
To further simplify the formulation, we now take the γ → 0 limit such that solvent density is no longer an independent component of the system but is instead determined by In other words, ∆ρ = 0 and, equivalently, The remaining system components {ρ p , ρ + , ρ − } are further constrained by electric charge neutrality of the system as a whole, i.e. σ p ρ p + ρ + − ρ − = 0, where σ p is the net electric charge of a polymer chain. To make this condition manifest, we introduce the overall salt density ρ 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, σ p ≥ 0, we have ρ + = ρ s and ρ − = σ p ρ p + ρ s ; for net-negative chains, σ p < 0, we have ρ − = ρ s and ρ + = −σ p ρ p + ρ s . Taking all the above considerations together and with k ≡ |k|, the MFT/RPA free energy per volume in units of accounts for translational entropy and follows from applying Stirling's approximation ln n! ≈ n ln n−n to the factorial prefactors in Eq. (19) and neglecting terms linear in ρ i arising from the nonlogarithmic −n part of the Stirling approximation because terms linear in ρ i in the free energy have no effect on phase separation properties. 38 In the second term on the right hand side of Eq. (34), φ p ≡ v b N ρ p /ρ 0 is the volume fraction occupied by polymer beads (and thus the volume fraction of the polymers themselves) and χ eff is a dimensionless effective Flory-Huggins χ-parameter that originates from the spatially short-range interactions [cf.
second term ofh in Eq. (29)]: Using Eq. (1), the summation over a in the above expression for χ eff can be reverted, i.e., 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-Williams 140 approximation used in Flory-Huggins theories of proteinlike heteropolymers. 137,141 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 ψ in Eq. (32) followed by subtraction of the free-energy contribution at ρ p = ρ s = 0, a subtraction that has no effect on phase separation properties. Given the free energy f (ρ p , ρ s ) in Eq. (34), phase diagrams may be constructed, for example, by matching the polymer and salt chemical potentials µ p,s = ∂f /∂ρ p,s and osmotic pressure Π = µ p ρ p + µ s ρ 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 χ 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 FTS because of its ability to afford, in principle, a full account of field fluctuations. 142 In FTS, each field (denoted generically as φ) 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 ∆x. The fields evolve in a fictitious complex-Langevin (CL) time t according to whereη φ (r, t) represents real-valued Gaussian noise with zero mean, i.e., Thermal averages of the real-space system can 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. 143,144 This approach has its origin in the development, beginning in the 1980s, of stochastic quantization 145 as a method for studying quantum field theories and their regularization [146][147][148] and has since been applied extensively to study properties of polymer solutions. 58 In the present model, the functional derivatives of the field Hamiltonian H in Eq.
whereρ i (r) in the first relation here in Eq. (39) is a field operator corresponding to the density of component i, and is given, respectively, bỹ for ions and solvents, whereas the polymer bead density operator,ρ b (r), is computed through forward and backward chain propagators q F,B (r, α) which are constructed by applying the following relations iteratively: 62,78,91 where Φ(r) = e −3r 2 /2b 2 /(2πb 2 /3) 3/2 , and initiating from q F (r, 1) = e −W 1 (r) and q B (r, N ) = e −W N (r) . Given q F,B , the polymer bead density operator can be computed as with In the second relation in Eq. (39),c(r) =c b (r) +ρ + (r) −ρ − (r) is a field operator corresponding to the electric charge density, wherein the polymer bead contributionc b (r) is given byc In the third relation in Eq. (39), the field operators for the nonelectric charge densities associated with the spatially short-range interactions,h a (r), are similarly given bỹ It can be readily verified that the position-independent solution to δH/δφ(r) = 0 for φ = η, ψ, ϕ a (where the expressions for the functional derivatives are given in Eq. (39)) yields exactly the MFT solution in Eq. (28) for position-independent η and ϕ a (together with any value of position-independent ψ): For any set of values of the position-independent fields, the field operators for number and charge densities become equal to their bulk counterparts, i.e.ρ b,±,w (r) = ρ b,±,w ,c(r) = ρ p σ p + ρ + − ρ − = 0 andh a (r) = ρ 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 h (r), where r = |r|, not only its spatial integral drV h (r) =V h (0) that appears as part of an effective Flory-Huggins χ-parameter, χ eff , in Eq. (36). Several functional forms for spatially short-range interactions have been used extensively in recent coarse-grained simulations of IDP LLPS. 40,43,117 These include the common Lennard-Jones (LJ) potential where E is the depth of the potential well and σ is a length scale, and LJ variations such as the shifted Weeks-Chandler-Andersen (WCA) potential 149 for finite-range purely repulsive interactions as well as the recently proposed Wang-Frenkel potential 150 with a finite spatial range r c , for r ≤ r c and V WF (µ, ν, r c , σ|r) = 0 for r > r c with 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 c . Note that the original notation in ref. 150 for the Wang-Frenkel form is largely followed in Eqs. (47) and (48). The symbols α WF and µ here should not be confused with the polymer bead label where interaction strength is represented by l h and the spatial range of the interaction is characterized by the reciprocal κ −1 of the screening coefficient κ (Fig. 4). The Fourier Pair Correlation Functions. Structural information of the molecular species in our model system can be inferred from pair correlation functions, 62,92 defined here as 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 models 62 and have been applied more recently to analyze experimental data on the biomolecular structure of nucleolar condensates. 153 The pair correlation functions are, in turn, related to the potentials of mean force (PMFs), U p,q (r), through the standard formula where the normalization factor ρ p ρ 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: In accordance with a common procedure in field theory, 154 these relations can be derived by adding source terms p drρ p (r)J p (r) to the particle HamiltonianĤ in Eq. (7) and then taking functional derivatives with respect to J p (r)s in the resulting field theory and finally setting all J p → 0. The formula for the self-correlation function ρ pρp may be obtained 62 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, 62,66,91 coarse-grained explicit-chain molecular dynamics (MD) simulations are conducted here to elucidate the ramifications of several proposed interaction schemes 40,66,111,112,116,117 and to assess the accuracy of our approximate analytical theories and FTS in capturing the essential physics of biomolecular LLPS. As before, 40,43 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, 91 for simplicity we will limit the present study to implicit-solvent simulations. 40,66 As the general features of the model can be found in the literature, 40,66,111,112,116,117 only an outline is provided here. Briefly, as in the above field-theory formulation, let n 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, . . . , n p and the residues (beads) along a polypeptide chain labeled by α or β = 1, 2, . . . , N . The total potential energy U T is a function of the residue positions {R i,α } that may be written as where U bond is the bond-length term for chain connectivity: with R i,α;j,β ≡ |R i,α − R j,β |. Parameter values for the MD potential U bond in Eq. (54) and for U el , U hh below are different across the three interaction schemes-Mpipi, 117 Urry, 111 and FB 112 -that we apply here to the Ddx4 IDRs. 66 To facilitate comparisons with prior coarsegrained MD simulation results in the literature, the parameters in these original references are largely used with the corresponding interaction schemes in the present simulations unless stated otherwise. For U bond , K bond = 19.1917 kcal mol −1 Å −2 , l = 3.81 Å for Mpipi; 117 K bond = 9.56 kcal mol −1 Å −2 , l = 3.82 Å for Urry; 111 and K bond = 2.39 kcal mol −1 Å −2 , l = 3.8 Å for FB. 112 The U el term in Eq. (53) is the electrostatic interaction: where, as in Eq. (13), σ i,α is the electric charge, in units of proton charge e, of the αth residue along the ith chain (as stated above for the field-theoretic formulation, σ i,α = σ α is independent of i in the present MD simulations), and κ D is the reciprocal of the Debye screening length. Mpipi 117 uses σ α = +0.75 for the positively charged residues arginine and lysine, σ α = +0.375 for histidine (to account for its incomplete ionization under physiological conditions), and σ α = −0.75 for the negatively charged residues aspartic and glutamic acids; Urry 111 uses σ α = +1 for arginine and lysine, σ α = 0 for histidine, and σ α = −1 for aspartic and glutamic acids; and FB 112 uses σ α = +1 for arginine and lysine, σ α = +0.5 for histidine, and σ α = −1 for aspartic and glutamic acids. All other residues have σ α = 0 in all of these three interaction schemes. In all of the present coarse-grained MD simulations, we set κ D = 0.1 Å −1 (i.e., a Debye screening length of 10 Å corresponding to a physiological monovalent salt concentration of ∼ 100 mM), apply a large-distance cutoff of 35 Å to the electrostatic interaction, and use a dielectric constant r = 40. Employing a dielectric constant moderately lower than the r ≈ 80 value of bulk water in coarse-grained MD simulation of IDP LLPS has recently been rationalized from a theoretical perspective 66,91 and this modeling practice is apparently not inconsistent with experimental observations of slower water orientational dynamics in living matter. 155 For the spatially short-range interaction term U hh in Eq. (53), Mpipi 117 uses the Wang-Frenkel form, 150 viz., where, by definition, V WF (E α,β , µ α,β , 1, (r c ) α,β , σ α,β |R i,α;j,β ) vanishes for R i,α;j,β > (r c ) α,β as specified in Eqs. (47) and (48). For all amino acid residue pairs α, β, Mpipi sets the distance cutoff (r c ) α,β = 3σ α,β and assigns ν = 1. Values of the interaction strength E α,β , length scale σ α,β , and the parameter µ α,β for the Wang-Frenkel functional form depend on the amino acid residue types r α and r β for residues α and β. The Mpipi values for E r,r , σ r,r , and µ r,r , where r and r each represents one of the twenty amino acid types, are provided in Supplementary where is in the Ashbaugh-Hatch form 156 with the function V LJ given by Eq. (46). For both Urry and FB, E α,β = E = 0.2 kcal mol −1 irrespective of the residue types r α , r β of residues α, β, and the length scale σ α,β = [σ(r α ) + σ(r β )]/2, where σ(r α ) is the diameter of the bead model for amino acid residue type r α given by the quantity σ in Table S1 of ref. 40. In both interaction schemes, λ hh α,β = [λ(r α )+λ(r β )]/2, where λ(r α ) is the hydrophobicity/hydropathy parameter for residue type r α . The λ(r α ) values for the Urry interaction scheme are given by the "Urry et al. normalized hydropathy scale" column of Table S2 in ref. 111, whereas the corresponding λ(r α ) values for the FB interaction scheme are provided in Table S7 of ref. 112. In simulations that use the Urry and FB interaction schemes, a cutoff distance of 20 Å is applied to the (U hh ) Urry/FB i,α;j,β 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, both with zero electric charge, but possess different hydrophobicities. The hydrophobic-polar patterns of these model sequences-commonly referred to as HP sequences in the protein literature 61,141,157,158 -are taken from recent studies of their LLPS. 132,133 In keeping with tradition of surfactant models, however, hydrophobic and polar beads were labeled, respectively, as "T" (tail) and "H" (head) in refs. 132, 133. 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 a test case, we use the KH interaction scheme 40 described in ref. 66 for the present LS-sequence simulations. Specifically, the spatially short-range interaction U hh is in the form of Eqs. (57) and (58)  MD Simulation of LLPS. Based on the models described above, we follow the recently developed, widely-applied slab method for simulation of IDP LLPS. 40,43,78 Because this simulation protocol, its rationale, and the procedure to construct phase diagrams from equilibrated simulations have been detailed elsewhere, 40,43,67 a brief outline here will suffice. As in our previous studies, 43,62,66,78,91 the GPU version of HOOMD-Blue software 159,160 is employed in the present simulations. For each simulation of the wildtype or a variant of the 236-residue Ddx4 IDR, 66 100 chains are randomly placed in a sufficiently large cubic box of dimension 300 × 300 × 300 Å 3 initially. Energy minimization is then performed to remove any steric clashes among the model molecules. This is followed by N P T 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 barostat 161,162 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. 163 After the initial N P T 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 × 150 × 150 Å 3 . At this point, the system is expanded along one of the axes (referred to as the z-axis) to 3000 Å for 10 ns with the temperature fixed at 100 K. Next, an N V T equilibration is performed for 2 µ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 µ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 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 × 130 × 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 based upon the success and limitation of the following examinations. We do so by using our formulation on four well-studied Ddx4 IDRs 9,19,66,98,111,117 as well as three selected hydrophobic-polar model sequences that have the same composition but different hydrophobic/polar sequence patterns. 132

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-tolysine (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 ε r,r s from the five interaction schemes illustrated in Fig. 2. In accordance with Eqs. (36) and (37) for χ eff , all nonzero eigenvalues λ a s are included in the application of the MFT/RPA Eq. (34) to the four Ddx4 IDRs. Fig. 5 shows the nonelectric charge profiles of the WT, RtoK and FtoA Ddx4 IDRs for the two most dominant eigenvalues λ 1,2 in each of the five interaction schemes considered. Significant differences in nonelectric charge profiles are observed across different interaction schemes (Fig. 5a- varible α) 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 may be 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 χ eff parameters, which may be computed from the eigenvalues λ a of the spectral decomposition of ε r,r in conjunction with the nonelectric charges q (a) rα using Eq. (36) or directly from the interaction matrix ε r,r itself via the identity in Eq. (37). In the present MFT/RPA formulation, which for simplicity assigns the same V h (r) to every residue-residue interaction scheme considered, the ratios between these χ eff parameters depend only on ε r,r but not V 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 (59) wherein WT and CS share the same χ eff because the WT and CS sequences have the same amino acid composition. Utilizing the ε r,r values of the interaction schemes 40,111,112,117 in Eqs. (36) and (37) leads to the following ratio: Phase diagrams can readily be constructed 78 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 byV h (0), is adjusted for each interaction scheme separately such that χ here, electric charge is +e for residues R and K, −e for D and E, and zero for all the other residues. Fig. 6 shows phase diagrams in the (ρ p , ρ 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 shows (ρ p , T ) phase diagrams of the same systems at a constant bulk protein density of ρ p = 10 mM (indicated by the vertical grey dashed lines), and at two bulk salt densities ρ 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 K ≤ T ≤ 510 K. At each temperature, we check if the point (ρ p , ρ s ) lies in the two-phase region and then identify the tie-line it intersects. The densities of the coexisting 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 ρ p = 10 mM. The coexisting protein densities (i.e., the phase dia-  Fig. 6) at multiple T . For a given T , the two coexisting ρ p s at bulk salt density ρ s = 100 or 300 mM are obtained as the ρ p values at the two ends of a (ρ p , ρ s ) coexistence tie-line that passes through the intersection point of ρ p = 10 mM (indicated by the vertical dotted lines in the present figure) and ρ s = 100 mM or ρ s = 300 mM, respectively, in the ρ p -ρ s plane. The color code for the Ddx4 variants are the same as that in Fig. 6. Dashed curves are used here for FtoA to facilitate visual discrimination between the near-overlapping coexisting curves for FtoA (red) and RtoK (green) in the Urry interaction scheme. gram) 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. 66,111,117 Two-solutespecies (ρ p , ρ 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 that of 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, 66,117 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. 9,19,98 For the RtoK and FtoA variants, the trend is underpinned by the effective Flory χ parameter ratios χ  (Table S7 for optimized FB-HPS parameters in ref. 112).
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 ∼ 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

Field-Theoretic Simulations of Phase Behaviors of Model Hydrophobic-Polar
Sequences with the Same Composition but Different Sequence Patterns. As a first 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 hydrophobicpolar sequences with the same 3/2 hydrophobic/polar composition but different sequence patterns of the hydrophobic and polar residues. We adopt three 20mer (N = 20) sequences with f T = 0.6 from Statt et al. 132 (hydrophobic composition is denoted as f T in this reference), namely their "[HT] 2 TH[T 2 H] 4 TH", "T 3 H 3 T 3 H 2 T 3 H 3 T 3 ", and "T 12 H 8 " sequences listed, respectively, as the last sequence in Fig. S6, the sequence in Fig. 5, and the "micelles" sequence in Fig. 6 in ref. 132. 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). For computational efficiency, we opted to use these short LS sequences without electrostatic interactions instead of the Ddx4 IDR sequences that are more than ten times longer for the first test of our FTS formulation. FTS of longer sequences with more complex interactions entails technical issues with regard to equilibration that remains to be tackled.
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 hydrophobicpolar (HP) potential 157 whereby only hydrophobic-hydrophobic interactions are favorable (governed by a LJ potential in the form of 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 potential 149 ). Within this HP interaction scheme, the LLPS critical temperatures (T cr = T c in their notation) for the first and second sequences are, in model units, T cr = 0.808 and T cr = 1.083, whereas the third sequence apparently does not have a critical point (Table S1 of ref. 132). Since we will be applying the KH interaction scheme to the LS sequences, with E rα,r β = −0.771, 0, and −0.304, respectively, for (r α , r β ) = (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. The present choice of using the KH potential for leucine and serine-with nonzero interactions for all three possible interaction pairs instead of the HP interaction scheme which is diagonalized by construction with only one nonzero interaction-is motivated by our desire to test the new FTS formulation with a more complex interaction scheme that requires diagonalization of the ε r,r matrix.
The FTS is performed on an M = 48 lattice with lattice spacing ∆x = a s = b/ √ 6 ≈ 1.6 Å, which amounts to a simulation box with a linear dimension of ≈ 75 Å. Using the KH interaction matrix, we truncate the summations over eigenvalues atã = 3 [Eqs. (1) and (2)]. The range of spatially short-range interactions governed by V h (r) in Eq. (49) is set to κ −1 = 5 Å and we let v b = v w = 1 for simplicity. All simulations are performed without salt, and since the LS sequences carry no electric charge, the field ψ(r) in Eq. (20) can be omitted. The total density of the system is set to ρ 0 = ρ w + N ρ p = 55.5 M, with a polymeric (LS sequence) component of ρ p = 300 mM giving a polymer bead volume fraction of φ = 0.11. The system thus contains n p ≈ 75 LS polymer chains and n w ≈ 1.2 × 10 4 solvent molecules. We set the compressibility to γ = 0.1b FTS of exactly incompressible systems. 91 However, we note that incompressible systems can alternatively be explored in FTS using the partial saddle-point approximation of ref. 164. Simulations are performed at l h = 0.5 Å (χ eff = 2.16) and at l h = 0.03 Å (χ 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 ∆t = 5 × 10 −4 b 3 using a semi-implicit integration scheme 165 generalized to multiple fields. 78 We perform 21 independent simulations for each sequence and parameter set, and use Q p,w to monitor system equilibration. At l h = 0.5 Å, the system requires 4×10 5 time-steps for equilibration, after which the field configuration are sampled every 50th step during the production run consisting of another 4×10 5 CL steps. At l h = 0.03 Å, 2×10 4 integration steps are sufficient for equilibration, followed by sampling every 50th step during a production run of 3 × 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 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. Colored curves in Fig. 8 correspond to l 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 p,p (r) ∼ 0 (Fig. 8b) at separations r larger than the condensate diameter, combined with polymer-solvent anticorrelation 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 e Wα(r) q F (r, α)q B (r, α)δ rα,X , X = L (leucine), S (serine), for the number density of L and S residues (the Kronecker δ r,r = 1 for r = r and δ r,r = 0 for r = r ). To facilitate visualization of the condensate structure, the snapshots in Fig. 9 are averaged over 10 3 consecutive CL time steps. Fig. 10 provides 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. 132. A similar behavior has been observed in Monte-Carlo simulations of 10mer hydrophobic-polar sequences, where the sequence of alternating hydrophobic-polar residues underwent LLPS whereas the diblock condensate exhibited micellar structure instead. 61 Assessing Analytical and FTS Results by Coarse-Grained Explicit-Chain MD. We now proceed to further evaluate the effectiveness of the above analytical and FTS approaches by comparing their predictions with the corresponding results from coarsegrained 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, on one hand, the 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 h (r) in Eq. (49) with the same l h and κ 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 σ α,β for different types of amino acid pairs r α , r β in Eq. (56). On the other hand, it is computationally much more efficient to account for explicit solvent in field-theoretic approaches than in MD. 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. In future investigations, it will be interesting to explore possible tuning of energetic parameters in our analytical and FTS formulations to optimize agreement with coarse-grained explicit-chain MD and/or experimental data to enable more broadly practical applications of field-theoretic techniques.
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, 117 Urry, 111 and FB 112 interaction schemes (Fig. 11). Corresponding coarse-grained explicit-chain MD phase diagrams using KH and HPS have been provided, respectively, in Fig. 4 and Fig. 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 sequencedependent LLPS of IDPs. In particular, in agreement with Fig. 6b of the original Mpipi study 117 (dotted curves in Fig. 11a, which were computed by using κ −1 D = 7.95 Å, r = 80 instead of the κ −1 D = 10 Å, r = 40 values used for the solid curves in Fig. 11a here) and consistent with the present analytical MFT/RPA results in Fig. 7 at lower model temperatures T 250 K, Fig. 11a exhibits the experimentally correct WT > CS > FtoA > RtoK rank ordering of LLPS propensities, 9,19,98 though the Mpipi-computed difference in LLPS propensity between WT and CS (∼ 6-10 K difference in T cr ) 117 is significantly smaller than the experimental estimation of ∼ 80 K difference in T cr at [NaCl] = 100 mM (ref. 19). For Urry, the LLPS propensity rank ordering WT > FtoA ≈ 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. A comparison of the global critical temperatures computed by analytical MFT/RPA and those simulated using coarse-grained MD at a fixed monovalent salt concentration of ∼ 100 mM for the interaction schemes Mpipi, Urry, and FB is provided in Fig. 12.
As discussed in previous works, 66,111,112,117 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 E R,K ≡ 20 r=1 E r,r /20 over the 20 amino acid types, where r is the label for arginine (R) or lysine (K). Entries for E r,r in Supplementary Table 11 of ref. 117 yield E R = 0.2212 kcal mol −1 and E K = 0.0563 kcal mol −1 , thus ∆E R−K ≡ E R − E 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 ∆(λ hh E) R−K = 0.0829 kcal mol −1 , which is significantly smaller. Here, for KH, ∆(λ hh E) R−K ≡ (λ hh E) R − (λ hh E) K , with λ hh E R,K ≡ − 20 r=1 λ hh r,r E r,r /20, and −λ hh r,r E r,r being the entries in Table S3 of ref. 40. For HPS, Urry, and FB, the corresponding difference in interaction strength is given by E∆λ R−K /2 ≡ E[λ(R) − λ(K)]/2 where E = 0.2 kcal mol −1 and λ(R) and λ(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  As noted above, the difference in LLPS propensity between WT and CS predicted by Mpipi is small, with ∆(T cr ) WT−CS ≈ 6 K in the original simulation 117 and ≈ 15 K in Fig. 11a, both much less than the experimentally estimated ∆(T cr ) WT−CS ∼ 80 K (refs. 19). In fact, the slightly higher ∆(T cr ) WT−CS we obtain is probably due in large part to the smaller dielectric constant r = 40 used for Fig. 11a instead of the r = 80 value used in the original study. By comparison, larger ∆(T cr ) WT−CS values are predicted by most of the other interaction schemes considered here: For HPS, ∆(T cr ) WT−CS ≈ 14 K for r = 80 (Fig. 3B  of ref. 66). For KH, ∆(T cr ) WT−CS ≈ 21 and 46 K, respectively, for r = 80 and 40 (Fig. 4  of ref. 66). For Urry and FB simulated here at r = 40, ∆(T cr ) WT−CS ≈ 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 histidine 117 instead of using the full proton charge e for electrostatic interactions among R, K, D, and E as in  Fig. 11 and the global T cr in MFT/RPA mostly agree for each of the three interaction schemes. Note that the discrepancies in the magnitudes of MD-and MFT/RPA-predicted T cr values shown here are largely a consequence of different energy scales used for MD and MFT/RPA. As stated in the discussion above for the results in Fig. 6, in order to facilitate comparison across different interaction schemes, the chosen energy scales for MFT/RPA are such that χ eff is always equal to 0.5 at a reference temperature T = 300K and l B = 7 Å for WT Ddx4 for all interaction schemes (thus all the blue lines for MFT/RPA in (b) are at the same level). In contrast, we largely follow the full interaction potentials provided in refs. 111, 112 and 117 for the Urry, FB, and Mpipi interaction schemes. It follows that only the rank orderings-but not the magnitudes-of T cr s for the Ddx4 variants are comparable in (b). 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 ±0.75e electric charges instead of ±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 κ −1 D = 7.95 Å instead of the κ −1 D = 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 ∆(T cr ) WT−CS it predicts. Secondary to electrostatics, another possible origin of the observed variation in ∆(T cr ) WT−CS predicted by different interaction schemes is their different spatially shortrange 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 ↔ R), 5 (D ↔ R), 2 (E ↔ K), and 1 (D ↔ 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 ∆(T cr ) WT−CS and this effect clearly depends on the interaction scheme. Now, the energetic effects of the (E ↔ R), (D ↔ R), (E ↔ K), and (D ↔ K) swaps may be characterized by E∆λ E−R /2, E∆λ D−R /2, E∆λ E−K /2, and E∆λ D−K /2; but these quantities can have different magnitudes and even different signs for different interaction schemes, as exemplified by the E∆λ E−R /2, E∆λ D−R /2, E∆λ E−K /2, and E∆λ 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. 13 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 LS1 observed in Fig. 13a 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, but the MD phase diagrams show that LS2 has a slightly higher LLPS propensity than LS1. Such minor differences are not unexpected 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.
The critical temperatures of the LS sequences shown in Fig. 13a 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. 13b-d). This observation is notable because LLPS propensities tend to increase with chain length 36 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. 13d 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. 132, 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. 132.
A sequence hydropathy decoration (SHD) parameter, defined as SHD = N −1 N α<β [λ(r α ) + λ(r β )]|α − β| −1 , was proposed recently as a predictor of prop- where ε r,r 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, SHD hh = 1.12, 1.05, and 1.17, respectively. As the SHD hh value of LS3 with higher LLPS propensity is larger than those of LS1 and LS2, this result suggests that SHD or SHD hh may be used as a predictor for LLPS as well. Further effort will be needed to examine whether the correlation between SHD or SHD hh and LLPS propensity of heteropolymers with spatially short-range interactions is as strong as that between sequence charge decoration (SCD) 86 and LLPS propensity of polyampholytes. 39,41,43

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 π-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, [166][167][168] future effort should aim to extend our formulation to other functional forms for spatially short-range potentials, 151 including temperature-dependent effects. 139,167 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.