A Simple Explicit-Solvent Model of Polyampholyte Phase Behaviors and its Ramifications for Dielectric Effects in Biomolecular Condensates
Jonas Wessén, Tanmoy Pal, Suman Das, Yi-Hsuan Lin, Hue Sun Chan
11February 6, 2021
A Simple Explicit-Solvent Model of PolyampholytePhase Behaviors and its Ramifications forDielectric Effects in Biomolecular Condensates
Jonas W
ESS´EN , , † Tanmoy P AL , , † Suman D AS , , † Yi-Hsuan L IN , , and Hue Sun C
HAN , ∗ Department of Biochemistry, University of Toronto, Toronto, Ontario M5S 1A8, Canada; Molecular Medicine, Hospital for Sick Children, Toronto, Ontario M5G 0A4, Canada † Contributed equally to this work. ∗ Corresponding authorE-mail: [email protected]; Tel: (416)978-2697; Fax: (416)978-8548Mailing address:Department of Biochemistry, University of Toronto, Medical Sciences Building – 5th Fl.,1 King’s College Circle, Toronto, Ontario M5S 1A8, Canada. a r X i v : . [ q - b i o . B M ] F e b Abstract
Biomolecular condensates such as membraneless organelles, underpinned by liquid-liquidphase separation (LLPS), are important for physiological function, with electrostatics—among other interaction types—being a prominent force in their assembly. Charge interac-tions of intrinsically disordered proteins (IDPs) and other biomolecules are sensitive to theaqueous dielectric environment. Because the relative permittivity of protein is significantlylower than that of water, the interior of an IDP condensate is a relatively low-dielectricregime, which, aside from its possible functional effects on client molecules, should facili-tate stronger electrostatic interactions among the scaffold IDPs. To gain insight into thisLLPS-induced dielectric heterogeneity, addressing in particular whether a low-dielectriccondensed phase entails more favorable LLPS than that posited by assuming IDP electro-static interactions are uniformly modulated by the higher dielectric constant of the puresolvent, we consider a simplified multiple-chain model of polyampholytes immersed in ex-plicit solvents that are either polarizable or possess a permanent dipole. Notably, simulatedphase behaviors of these systems exhibit only minor to moderate differences from thoseobtained using implicit-solvent models with a uniform relative permittivity equals to thatof pure solvent. Buttressed by theoretical treatments developed here using random phaseapproximation and polymer field-theoretic simulations, these observations indicate a partialcompensation of effects between favorable solvent-mediated interactions among the polyam-pholytes in the condensed phase and favorable polyampholyte-solvent interactions in thedilute phase, often netting only a minor enhancement of overall LLPS propensity from thevery dielectric heterogeneity that arises from the LLPS itself. Further ramifications of thisprinciple are discussed.
INTRODUCTION
Biomolecular condensates, —physico-chemically based in large measure upon liquid-liquid phase separation (LLPS) underpinned by multivalent interactions among intrinsi-cally disordered proteins (IDPs), proteins containing intrinsically disordered regions (IDRs),folded domains, folded proteins, and nucleic acids —are now widely recognized as a meansfor living organisms, eukaryotes as well as prokaryotes, to functionally compartmentalizetheir internal space. Among their multifaceted physiological roles, the droplet-likecondensates, especially the intracellular varieties referred to as membraneless organelles,provide specialized, tunable micro-environments, distinct from their surroundings, that fa-cilitate certain biochemical processes.
One of the many possible mechanisms by whichthey serve this task is by preferential partitioning of certain reactants into the condensates.Examples include retention of selected types of small nuclear ribonucleoproteins in Cajalbodies, a function likely plays a critical role in pre-mRNA splicing, partitioning of single-stranded but not double-stranded DNA into condensates of RNA helicase Ddx4 that mimicnuage membraneless organelles—thus serving probably as a molecular filter for functionalnucleic acid processing, and concentration of transcription regulators in TAZ (tafazzin)nuclear condensates to promote gene expression. LLPS of IDPs and IDRs is known to underlie an increasing number of biomolecularcondensates, as exemplified by well-studied condensates of N-terminal IDR domains of,respectively, the DEAD-box RNA helicase 4 (Ddx4), the RNA-binding protein fusedin sacroma (FUS), and the P-granule DEAD-box RNA helicase LAF-1. Comparedto amino acid sequences encoding for globular proteins, IDPs and IDRs are depleted ofhydrophobic but enriched with polar and charged residues. Consistent with these trends,electrostatics is one of the major driving forces for biomolecular LLPS, together withcation- π and π - π interactions which can be equally or sometimes even more important, and hydrophobic interactions which can be important in some cases, particularly forsystems that exhibit a lower critical solution temperature rather than an upper criticalsolution temperature. As far as electrostatics is concerned, experiments on a chargescrambled mutant of Ddx4, phosphorylated variants of FUS, phosphorylated variants ofan IDR of the fragile X mental retardation protein (FMRP), and charge shuffled variantsof LAF-1 indicate that effects of charged interactions on LLPS properties depend not onlyon an IDP/IDR’s net charge but also its sequence charge patterns. Recently, some ofthese and related sequence-dependent effects on LLPS have been modeled successfully andrationalized by polymer theory using random phase approximation (RPA), restrictedprimitive models, Kuhn length renormalization, coarse-grained explicit-chain moleculardynamics (MD) and Monte Carlo simulations, and field-theoretic simulations. In these recent theoretical and computational studies of electrostatic effects in LLPS, therelative permittivity, (cid:15) r , for the solvent-mediated intra- and interchain electrostatic inter-actions among charged IDPs are often taken to be uniform and equal to that of bulk water( (cid:15) r ≈ or left unspecified as a constant (cid:15) r that is subsumed in the Bjer-rum length or a reduced temperature, or as a fitting parameter for matchingtheory to experiment, yielding constant (cid:15) r values in the range of ≈ A commonfeature of these formulations is that while the value of (cid:15) r may differ from model to model, (cid:15) r is always assumed to be a position-independent constant for a given model irrespectiveof the local IDP density. However, proteins and nucleic acids have substantially lowerrelative permittivities than bulk water. Inside the droplet, the density of water is reducedbecause of the higher concentrations of protein and nucleic acids. Moreover, the watermolecules that remain, being likely near protein surfaces, generally entail a smaller rela-tive permittivity than that of bulk water due to restrictions on their rotational freedom. Taken together, these considerations lead one to expect that the effective permittivity inside a condensed-phase IDP droplet is lower than that of the dilute phase. Such a dif-ferent dielectric environment could affect the behaviors of proteins and nucleic acids as well as the partitioning and sequestration of a variety of “client” molecules in the con-densed droplet and thus contribute to the maintenance of biomolecular condensates andthe functional biophysical and biochemical processes within them. The significantly different dielectric environments expected for the dilute phase on onehand and the condensed phase on the other raise basic questions regarding the aforemen-tioned common practice of using a single, constant (cid:15) r —rather than a variable (cid:15) r that dependson local IDP density—for the electrostatic interactions in LLPS theories and computationalmodels. One obvious concern is that because the effective (cid:15) r inside the condensed phase issubstantially lower than that of water, the LLPS-driving electrostatic interactions in thecondensed phase should be modulated by this lower value of (cid:15) r and therefore much strongerthan that posited by assuming a uniform (cid:15) r of bulk water. This consideration suggeststhat the propensity to phase separate can be much higher than that predicted by assuminga constant (cid:15) r . In apparent agreement with this intuition, earlier studies using modifiedRPA formulations with an (cid:15) r ( φ ) that depends on IDP volume fraction φ show significantincreases in LLPS propensities. Subsequently, however, a more nuanced analysis thattakes into account a “self-energy” term in the (cid:15) r → (cid:15) r ( φ ) RPA formulation indicates thatthere is a trade-off between the enhanced effective electrostatic interactions among the IDPsin the condensed phase (which tend to increase LLPS propensity) and the enhanced inter-actions between the IDPs and the solvent in the dilute phase (which tend to decrease LLPSpropensity), resulting in a relative minor overall effect on LLPS propensity. While the insights provided by these preliminary investigations are valuable—not theleast for formulating the questions needed to be tackled, our understanding of the effect ofthe biomolecule-dependent dielectric environment on the energetics of LLPS remains quitelimited. This is because the (cid:15) r ( φ ) in these formulations are obtained from effective mediumtheories based on an implicit treatment of water rather than having dielectric prop-erties emerge from an explicit model of water. Indeed, a shortcoming of this approach isthat predictions based on different effective medium theories can differ substantially. With this in mind, the present work aims to attain a more definitive physical picture bymodeling sequence-dependent LLPS of polyampholytes with an explicit polar solvent. Re-cently, new insights have been gained into the interactions between salt ions and proteinas well as the diffusive dynamics in condensed FUS and LAF-1 IDR droplets by usingthe TIP4P/2005 model of water. The explicit water model, however, is absent in theinitial thermodynamic simulation to assemble the condensed phase. Separately, explicitcoarse-grained water has been applied successfully to simulate condensates of FUS IDR using the MARTINI coarse-grained force field, though polar versions of the coarse-grainedwater model is not employed. Since our goal is to better understand the general phys-ical principles rather than to achieve quantitative match with experiment, here we takea complementary approach by using simple dipoles—each consisting only of one positiveand one negative charge—as a model for a water-like polar solvent. In this regard, ourconstruct shares similarities with a simple lattice dipole model of water restricted to an es-sentially one-dimensional configuration as well as three-dimensional simple dipole modelsfor polar liquids. Unlike our previous implicit-solvent formulations, here the differentdielectric environments for the dilute and condensed phases arise physically from the polarnature of the dipoles rather than as a prescribed function of local solvent and IDP con-centrations alone. Importantly, for our purpose, the model’s simplicity allows equilibratedLLPS systems comprising of both explicit dipole solvent molecules and polyampholytes tobe simulated.Consistent with the trend observed in a recent implicit-solvent RPA analysis thatinclude a self-energy term, polyampholyte phase properties simulated in explicit-chainmodels with the present explicit dipole model of polar solvent are quantitatively similarto—albeit not identical with—those simulated without the explicit solvent but with auniform relative permittivity corresponding to that of the pure polar solvent. Whether theLLPS propensity of the explicit-solvent system is slightly higher or slightly lower than thatof the corresponding implicit-solvent system likely depends on the system’s structural andenergetic details. Nonetheless, our results suggest, assuringly, that the error in employingan IDP-concentration-independent bulk-water (cid:15) r ≈
80 in LLPS simulations may be relatively small even when the dielectric environments inside and outside of thecondensed phase are significantly different. To address the generality of this putativeprinciple, we have also developed corresponding field-theory formulations for polyam-pholytes with a polar solvent modeled as permanent or polarizable dipoles that are treatedexplicitly in the theory. Similar permanent and polarizable dipole models wereused to study dielectric properties of polar solutions and the effect of electric field onphase stability of charged polymers; but have not been applied to tackle the questionat hand. Here we find that the trend of phase behaviors predicted using RPA of ourfield-theory formulation regarding explicit versus implicit solvent are in agreement withthat obtained by explicit-chain simulations, suggesting that a compensation between morefavorable polyampholyte-polar-solvent electrostatic interactions in a high- (cid:15) r dilute phaseand more favorable effective inter-polyampholyte electrostatic interactions in a low- (cid:15) r condensed phase is a rather robust feature of polyampholyte LLPS systems. Details of ourmethodology and findings are described below. MODELS AND METHODS
The formulation and computational methodology of our coarse-grained explicit-chainmodel of polyampholyte LLPS in the presence of explicit dipoles as a model polar solventis adapted from—and follows largely—a protocol developed recently for implicit-solventsimulations of LLPS of chain molecules, an approach our group has applied previouslyto study LLPS of model polyampholytes and to address experimental data on IDR LLPS. Following refs. 53, 55, but with a notation more in line with that in ref. 51, let n p be thetotal number of polyampholyte chains in the system which are labeled by µ, ν = 1 , , . . . , n p , N be the number of monomers (beads) in each chain, labeled by i, j = 1 , , . . . , N , and n w be the total number of dipoles (which may be viewed as serving as model water moleculesand hence the subscript “w”) which are labeled by υ = 1 , , . . . , n w . Note that a dipole inthe present formulation is equivalent to a polyampholytic dimer chain carrying a pair ofopposite charges, denoted here as q d = ±| q d | , on its two monomer beads.The total potential energy U T , expressed in the equation below as a sum of contributions,is a function of the positions of the monomers in the polyampholyte chains and in the dipoles: U T = U bond + U LJ + U el (1)with U bond being the bond-length term for chain connectivity that applies to the polyam-pholytes as well as the dipoles: U bond = K bond (cid:34) n p (cid:88) µ =1 N − (cid:88) i =1 ( r µi,µi +1 − a ) + n w (cid:88) υ =1 ( r υ + ,υ − − a ) (cid:35) (2)where, in general, r µi,νj is the distance between monomer i in polyampholyte chain µ andmonomer j in polyampholyte chain ν , r υ ± ,υ (cid:48) ± is the distance between the position of thepositive (+) or negative ( − ) charge on dipole υ and the position of the positive (+) ornegative ( − ) charge on dipole υ (cid:48) , and a is the reference bead-bead bond length that servesas the length scale of the model. Here we set K bond = 75 , ε/a , where ε is an energy scaleto be defined below. This value of K bond , which is in line with the TraPPE force field, entails a very stiff spring constant; thus the dipoles are essentially permanent (fixed)dipoles (with magnitude of dipole moment µ ≈ | q d | a ) because they are practically notpolarizable by fluctuation of the length ( ≈ a ) of the bond connecting a dipole’s + | q d | and −| q d | beads. q d is given in units of the elementary electronic charge, e , below. (Thesymbol µ for dipole moment magnitude is not to be confused with the polyampholyte label µ = 1 , , . . . , n p ). U LJ in the above Eq. 1 accounts for the non-bonded Lennard-Jones-like interactions,which in the present model takes different forms depending on whether the monomer beadsinvolved are both from polyampholyte chains (pp), both from dipoles (dd), or one from apolyampholyte chain and the other from a dipole (pd): U LJ = U LJ , pp + U LJ , dd + U LJ , pd . (3)As before, the polyampholyte-polyamopholyte term U LJ , pp is the standard 12-6 LJ poten-tial: U LJ , pp = 4 ε LJ n p (cid:88) µ,ν =1 N (cid:88) i,j =1( µ,i ) (cid:54) =( ν,j ) (cid:34)(cid:18) ar µi,νj (cid:19) − (cid:18) ar µi,νj (cid:19) (cid:35) , (4)where ε LJ is the well depth and a is the bead diameter which is equal to the reference bondlength in Eq. 2 above. During simulation, this potential is truncated at r µi,νj = 6 a forcomputational efficiency (individual terms in U LJ , pp set to zero for r µi,νj ≥ a ). To avoidfreezing of the dipole solvent at relatively low temperatures and other complications as wellas to facilitate comparison with analytical theory, the dipole-dipole and polyampholyte-dipole contributions in Eq. 3 are chosen to be consisting of only excluded-volume (ex)repulsion terms, viz., U LJ , dd = n w − (cid:88) υ =1 n w (cid:88) υ (cid:48) = υ +1 [ ˜ U exLJ ( r υ + ,υ (cid:48) + ) + ˜ U exLJ ( r υ − ,υ (cid:48) − ) + ˜ U exLJ ( r υ + ,υ (cid:48) − ) + ˜ U exLJ ( r υ − ,υ (cid:48) + )] (5)and U LJ , pd = n p (cid:88) µ =1 N (cid:88) i =1 n w (cid:88) υ =1 [ ˜ U exLJ ( r µi,υ + ) + ˜ U exLJ ( r µi,υ − )] , (6)where the individual purely repulsive terms in the above summations take the Weeks-Chandler-Andersen form: ˜ U exLJ ( r ) ≡ (cid:40) ˜ U LJ ( r ) + ε LJ , for r ≤ / a , for r > / a . (7)which is obtained, as specified here, by performing a cutoff and a shift on an individualstandard LJ term ˜ U LJ ( r ) ≡ ε LJ [( a/r ) − ( a/r ) ] in Eq. 4.In a similar vein as Eq. 3, the electrostatic interaction U el in Eq. 1 is written as a sumof pp, dd, and pd components: U el = U el , pp + U el , dd + U el , pd . (8)The polyampholyte-polyampholyte contribution is given by U el , pp = n p (cid:88) µ,ν =1 N (cid:88) i,j =1( µ,i ) (cid:54) =( ν,j ) σ i σ j e π(cid:15) (cid:15) r r µi,νj , (9)where σ i is the charge of the i th monomer bead in units of elementary electronic charge e ,( σ i depends on the polyampholyte sequence but is independent of µ, ν ), and (cid:15) is vacuumpermittivity. For simulation with explicit dipoles, relative permittivity (cid:15) r is set to unity forthe U el , pp term in Eq. 9; but (cid:15) r > U el , dd = n w − (cid:88) υ =1 n w (cid:88) υ (cid:48) = υ +1 q e π(cid:15) (cid:32) r υ + ,υ (cid:48) + + 1 r υ − ,υ (cid:48) − − r υ + ,υ (cid:48) − − r υ − ,υ (cid:48) + (cid:33) , (10)where the charges q d on the dipoles are in units of elementary electronic charge e , and U el , pd = n p (cid:88) µ =1 N (cid:88) i =1 n w (cid:88) υ =1 σ i q d e π(cid:15) (cid:32) r µi,υ + − r µi,υ − (cid:33) . (11)As in ref. 55, the energy scale ε introduced above in conjunction with K bond is chosen tobe the strength of electrostatic interactions at a monomer-monomer separation a for thepolyampholytes, i.e., ε = e / (4 π(cid:15) (cid:15) r a ), and the well depth ε LJ of LJ-like interactions in Eq. 4is set to be equal to ε , i.e., ε LJ = ε , as in Fig. 4a of ref. 55. These specifications completethe description of the total potential energy function U T in Eq. 1 used in our simulations.Fig. 1 depicts the N = 50 model polyampholyte sequences and the model dipole solventmolecule used in this study. The three overall neutral sequences (each containing 25 positiveand 25 negative beads), sv15, sv20, and sv30, were introduced by Das and Pappu. Thesesequences have been used in recent explicit-chain simulations of sequence-dependent IDPbinding and LLPS.
Here, they serve to assess whether the sequence-dependent LLPS ob-served in explicit-chain, implicit-solvent simulations holds for explicit-chain-simulatedLLPS in the presence of explicit dipole solvent molecules, bearing in mind that these se-quences cover a considerable variation of charge patterns as quantitatively characterizedby the intuitive blockiness parameter κ (ref. 44; κ = 0 . . .
0, respectively,for sv15, sv20, and sv30) and the theory-derived “sequence charge decoration” parameter (SCD; − SCD = 4 .
35, 7 .
37, and 27 .
8, respectively, for sv15, sv20, and sv30).
EEEEEEEEEEEEEEEEEEEEEEEE E KKKKKKKKKKKKKKKKKKKKKKKKKE EEEE E EEEE E EEEE E EEEE E EEEE KKKKK KKKKKKKKKK KKKKKKKKKKK KK EEEEEEEEEEEEEEEEEEEEEEEEE KKKKKKKKKKKKKKKKKKKKKK K sv15sv20sv30 (a) (b) AB dipole Fig. 1:
Model polyampholytes and dipole solvent molecules. (a) The N = 50 model sequences studiedin the present work. As in ref. 44, positively and negatively charged monomer beads are labeled as “K”(lysine) and “E” (glutamic acid), respectively, and depicted in blue and red. (b) Dipole as model polarsolvent molecule. Underscoring that the charges on the solvent dipoles can be different from those on themodel polyampholytes in (a), i.e., | q d | does not necessarily equal unity, the positively charged (blue, “B”)and negatively charged (red, “A”) beads in (b) do not share the “K” and “E” labels in (a). As in our previous studies, molecular dynamics simulations of the present modelare performed by the GPU-version of the HOOMD-blue simulation package using arecently developed methodology for simulating phase behaviors as well as binding ofIDPs. The computational protocol has been described in detail in previous works. Electrostatic interactions are treated using the PPPM method implemented in theHOOMD package. For our simulation of polymers (polyampholytes) and dipoles, n p = 100polymer chains are first randomly placed in a relatively large simulation box of dimensions40 a × a × a . This system is energy-minimized for a period of 500 τ with a timestep of0 . τ , where τ ≡ (cid:112) ma /ε with m being the mass of each monomer bead. For simplicity,similar to the models in ref. 55, all monomer beads making up the polyampholyte chains andthe dipole dimers are assigned the same mass. The initial equilibration is then performedfor a period of 1 , τ at T ∗ = 4 .
0. The reduced temperature T ∗ ≡ k B T /ε , where k B isBoltzmann constant and T is absolute temperature, contains a factor of (cid:15) r via the abovedefinition for ε . For direct comparisons of temperatures in different dielectric environmentson the same footing, however, it is often useful to use the vacuum reduced temperature ( T ∗ for (cid:15) r = 1), denoted in the present work as T ∗ ≡ T ∗ (cid:15) r = 4 π(cid:15) k B T ae . (12)0The velocity-Verlet algorithm is applied to propagate the equations of motion and periodicboundary conditions are applied to all three dimensions. The system is then compressedisotropically at T ∗ = 4 . , τ using linear scaling to a sufficiently higher density of ∼ . m/a by confining the polymers to a box of size 19 a × a × a . We next expand thebox 15 times along one of the spatial dimensions (labeled as z ) at a much lower temperatureof T ∗ = 0 .
5, resulting in final box dimensions of 19 a × a × a . After this procedure,the system is equilibrated for 1 , τ at the same T ∗ = 0 . . m/τ (ref. 95). The Packmol package is subsequently usedto insert n w = 12 ,
000 dipoles into the simulation box while keeping the positions of allthe polymer beads fixed. During the insertion process, a suitable distance criterion is setto avoid steric repulsion. After inserting the dipoles, we equilibrate the system again fora longer period of 30 , τ at several desired temperatures T ∗ ≥ . T ∗ < . T ∗ ≥ . , τ , during which snapshots aresaved every 10 τ for the determination of various ensemble properties and density profilesof the polymers and dipoles in the simulation box. Coexistence curves (phase diagrams) ofthe polymers are determined from the polymer density distributions as described before. Density of dipoles ( ρ w ) and density of polyampholytes ( ρ ) are reported by the numberdensities of monomer beads in the two types of molecules, respectively, in units of 1 /a .For the study of pure polymers, i.e., implicit-solvent systems of polyampholytes withno dipoles, the procedure is identical to that in ref. 55, viz., n p = 500 polymer chains areinitially placed randomly in a cubic box of length 70 a ; after energy minimization and initialequilibration, the 70 a -cubic box is compressed to a cubic box of length 33 a . The simulationbox for the system is then expanded eight times along one of the axes (named z as before) to reach a final dimensions of 33 a × a × a . As for the polymer plus dipole simulationsdescribed above, equilibration is performed, under periodic boundary conditions, in theexpanded box for 30 , τ , which is then followed by a production run for 100 , τ withsnapshots saved every 10 τ for subsequent analysis. Further details regarding simulationtemperatures and boundary conditions are provided in ref. 55.For the study of pure dipole solvent, i.e., systems that contain dipoles but no polymers,the main purpose of the simulations is to determine the effective relative permittivity,under various conditions, of bulk solvent as modeled by the dipoles. We use n w = 10 , ρ w ). We first focus on ρ w = 0 .
205 and ρ w = 0 . a and 43 a . The simulations are conducted for a broadrange of temperatures, namely T ∗ = 3 .
0, 4 .
0, 5 .
0, 6 .
0, and 7 .
0, as well as for a broad rangeof dipole moments, with q d = ± ± ± ±
4, and ±
5. To probe the dependence of theeffective relative permittivity on both density and temperature, we further fix q d = ± ρ w from 0 .
205 to ≈ .
55 at intervals of ≈ .
05 ( ρ w = 0 . . . . . . . . T ∗ = 3 . . ρ w as wellas a temperature-independent effective relative permittivity, we fix ρ w = 0 .
252 and usetheoretically fitted values of q d (see below) to conduct simulations at T ∗ = 3 .
0, 4 .
0, 5 . .
0, 7 .
0, and 8 . a ,are: 46 .
0, 43 .
0, 40 .
5, 38 .
5, 36 .
8, 35 .
4, 34 .
2, and 33 . τ , then equilibrate the system at the desired temperature using Langevin dynamicswith the same low friction coefficient of 0 . m/τ for 500 τ . The subsequent production runis performed for 4 , τ and snapshots are collected every 1 . τ for analysis. RESULTS
Effective Permittivities Of Pure Dipole Solvents.
We begin by determining base-line dielectric constants of the bulk solvents in our model by conducting simulations of puredipole systems under a variety of conditions as outlined above and utilizing the followingexpression for the effective relative permittivity: (cid:15) r = 1 + (cid:104) M T · M T (cid:105) − (cid:104) M T (cid:105) · (cid:104) M T (cid:105) (cid:15) k B T (13)where M T is the total dipole moment vector and Ω is the volume of the system, and (cid:104) . . . (cid:105) represents averaging over the simulated ensemble, effectuated here by averaging oversnapshots. The total dipole moment is given by M T = e (cid:80) i q i R i , where the summation isover monomer beads of the system/subsystem of interest, and q i and R i are, respectively,the charge and position vector of the i th bead. Our primary focus here is the dielectriccontribution from the dipole solvents, in which case the summation is over the beads ofthe dipole dimers. In situations where the contribution from the polyampholytes to thedielectric environment is also taken into consideration, the summation would encompass allmonomer beads in the system. (Note that the magnitude M T in Eq. S12 of ref. 34 for (cid:15) r should be replaced by the vector M T .)2 T *0 r ± ± ± ± ± (a) w = 0.205 T *0 ± ± ± ± ± (b) w = 0.252 Fig. 2:
Effective relative permittivity of pure dipoles as a model of polar solvent. Effective (cid:15) r is calculatedby Eq. 13 with ρ w = 0 .
205 (a) and 0 .
252 (b) for dipoles with the q d values indicated (shown as ±| q d | onthe right of the curves) as functions of reduced temperature T ∗ . Solid lines connecting the simulated datapoints (circles) are merely a guide for the eye. Dashed curves representing Eq. 14 for idealized uncorrelateddipoles are included for reference. Dotted curves show the empirical fit for the different | q d | values inaccordance with Eq. 15. Fig. 2 shows the variation of effective relative permittivity of pure dipole solvents withrespect to dipole density ( (cid:15) r for ρ w = 0 .
205 and 0 .
252 are compared), temperature T ∗ , andthe charges q d on the dipoles that make up the model polar solvent. As expected, since (cid:15) r ispositively correlated with the ease at which the solvent can be polarized, the simulated datain Fig. 2 indicate that (cid:15) r increases with increasing | q d | and ρ w . This is because the lattertwo increases correspond to having stronger charges on the dipoles (hence a larger dipolemoment µ ≈ | q d | a ) and the presence, in a given volume, of a larger number of dipoleswhose orientations can be perturbed by the insertion of other charges into the system.As temperature T ∗ increases, however, (cid:15) r decreases because biased distributions of dipoleorientations are disfavored by thermal agitation in general, as reflected, to a degree, by the1 /T dependence of the second term in Eq. 13.Inasmuch as the orientations of dipoles can be assumed to be uncorrelated, the analyticalformula (cid:15) r = 1 + 4 πl B µ (cid:16) ρ w (cid:17) , (14)where l B = e / (4 π(cid:15) k B T ) is the vacuum Bjerrum length and ( ρ w /
2) is the number densityof dipoles (because ρ w is defined to count two monomer beads per dipole as stated above),3applies for a collection of fixed dipoles with µ being the magnitude of the individual dipolemoment (see, e.g., refs. 89, 105). This formula does not provide a sufficiently accuratedescription of our simulation data (dashed curves in Fig. 2), however, because although themagnitude of the dipole moment for a given | q d | in our model may be considered practicallyfixed at µ = | q d | a because of the stiff spring constant K bond as noted above, the orientationsof the dipoles are correlated to a degree because of inter-dipole electrostatic interactions. Asexpected, the correlation among dipole orientations increases with increasing | q d | , leadingto larger mismatches seen in Fig. 2 between the simulated (cid:15) r data points and the dashedcurves representing Eq. 14. Nonetheless, we find that the following empirical formula, whichtakes a form similar to Eq. 14 but with an extra multiplicative factor linear in µ , viz., (cid:15) r = 1 + 4 πl B µ (cid:16) ρ w (cid:17) ( C + C µ ) (15)is capable of providing a reasonably good global fit for all the simulated (cid:15) r values in Fig. 2(dotted curves) when C = 1 . C = 0 . (cid:15) r values of dipole systems beyond those we have directly simulated. (a) (b) Fig. 3:
Dependence of effective relative permittivity of model dipole solvent on density, dipole moment,and temperature. (a) Variation of effective (cid:15) r (Eq. 13) with respect to dipole density ρ w and reducedtemperature T ∗ (marked on the right of the curves) for dipoles sharing the same | q d | = 3 .
0. The verticaldashed arrow across the curves indicates increasing temperature. (b) At constant ρ w = 0 . (cid:15) r (Eq. 13) for q d s estimated by Eq.15 ( q d values indicated beside simulated data points) to resultin an essentially constant (cid:15) r ≈ . T ∗ ). Note that the vertical scale covers only a narrow range of (cid:15) r from ≈ .
85 to ≈ .
10. Solidlines connecting the simulated data points (circles in (a) and (b)) are merely a guide for the eye.
Further data in Fig. 3a on the dependence of simulated (cid:15) r value on temperature anddipole density ρ w suggest an approximate linear relationship between (cid:15) r and ρ w . Consistent4with the observation mentioned above that the dipoles are well dispersed, and thereforefluid-like, for T ∗ ≥ .
0, the variation of (cid:15) r with respect to ρ w and T ∗ is seen to be smooth inthis set of simulations. Making use of this observation in Fig. 3a and the empirical Eq. 15,we obtain a set of µ values for dipole solvents that would maintain an essentially constant (cid:15) r over a range of different temperatures (Fig. 3b). Because (cid:15) r for dipoles with a constant µ decreases with increasing T ∗ (Fig. 2), the µ value has to increase with increasing T ∗ tomaintain an essentially temperature-independent (cid:15) r (Fig. 3b). We will use these µ valuesto construct temperature-independent solvent dielectric environments for the simulationstudy of polyampholyte LLPS below. Phase Behaviors Of Polyampholytes In Dipole Solvent.
With the baseline ef-fective relative permittivities of dipole models of polar solvent established, we now turn toLLPS behaviors of the polyampholyte chains immersed in these model solvent molecules.For most of these studies, we set the effective bulk-solvent (cid:15) r ≈ .
0. This relatively small ef-fective (cid:15) r value is chosen for computational tractability. We refrained from assigning highlypositive and negative charges to the model solvent molecules as well as having very highdensities of these molecules as the model solvent, both of these features would increaseeffective (cid:15) r but would decrease sampling efficiency. As it will become apparent below,the present model is useful for gaining insights into general physical principles regardingelectrostatics-driven LLPS in polar solvent environments, (cid:15) r ≈ . (cid:15) r ≈
80 for water notwithstanding.An overview of the properties of the present polyampholytes-and-dipoles systems is pro-vided in Fig. 4, which shows the density profiles of the polyampholytes and the dipoles forthe three sequences sv15, sv20, and sv30, each at four different simulation temperatures. Acondensed polyampholyte droplet, corresponding to a local peak in polymer bead density ρ (red curves), is seen in every panel of this figure. Expectedly, each of these droplets coincideswith a local minimum in dipole density (blue curves). As temperature increases from topto bottom in Fig. 4, effects of configurational entropy gain in significance. Consequently,the packing of polyampholytes in the droplet loosens up (lower and broadened peak, redcurves) and a concomitant relative increase in dipole density ensues (shallower but widerdip, blue curves). In other words, more solvent molecules reside in the polyampholyte-richphase as temperature increases.Fig. 4 shows further that the effective (cid:15) r contributed by the dipole solvent decreases insidethe polyampholyte-rich droplet (dips along green curves at z -positions of the droplets),tracking the decrease in ρ w . The z -dependent relative permittivity in Fig. 4 is computedusing the formulation in ref. 106 for “slab geometry” in which one spatial dimension ( z )is distinct from the other two spatial dimensions ( x, y ). The (cid:15) r ( z ) plotted in Fig. 4 is thecomponent of the diagonal permittivity tensor parallel to the latter two dimensions i.e., the5 x – x , or equivalently the y – y compoent, given by the expression (cid:15) r (cid:107) ( z ) = 1 + 12 (cid:15) k B T (cid:104) (cid:104) m (cid:107) ( z ) · M T (cid:107) (cid:105) − (cid:104) m (cid:107) ( z ) (cid:105) · (cid:104) M T (cid:107) (cid:105) (cid:105) (16)where the two-dimensional vectors m (cid:107) ( z ) and M T (cid:107) are both confined to the x – y plane. m (cid:107) ( z ) is the x – y component of the local dipole moment density (dipole moment within thelocal volume sampled divided by the local volume) and M T (cid:107) is the x – y component of thedipole moment of the entire system, and (cid:104) . . . (cid:105) denotes averaging over snapshots (Eq.13 ofref. 106). Note that this and related expressions for local (cid:15) r involves not only local but alsoglobal dipole moments, and Eq. 13 can be recovered from Eq. 16 by recognizing that m (cid:107) ( z ) → M T (cid:107) / Ω when the sampling volume for m (cid:107) ( z ) encompasses the total volume Ωand that when the system is isotropic after subtracting an overall average dipole moment (cid:104) M T (cid:107) (cid:105) , the relation (cid:104)| M T (cid:107) − (cid:104) M T (cid:107) (cid:105)| (cid:105) = (2 / (cid:104)| M T − (cid:104) M T (cid:105)| (cid:105) holds. In our calculationof m (cid:107) ( z ), we consider the dipoles in slices of thickness 10 a in the z -direction, and includeonly those dipoles whose centers of mass are inside the bin. The z – z component of thepermittivity tensor—the orthogonal component—may also be obtained and is expectedto yield values similar to that given by Eq. 16; but the more limited statistics availableto this calculation due to the restriction to one spatial dimension would likely presentcomputational challenges. For this reason we do not pursue it here, while recognizingthat Fig. 4, as it stands, illustrates reasonably well that the solvent dielectric environmentexperienced by polyampholyte chains inside the condensed droplet can be quite differentfrom that of the outside in the dilute phase.At temperatures significantly higher than those simulated in Fig. 4, condensed dropletsdisassemble. An example is provided by the snaphots in Fig. 5 for the polyampholytesequence sv20 at two temperatures. Consistent with the ρ w profiles in Fig. 4, at a moderatetemperature of T ∗ = 5 .
5, Fig. 5a–c show that although a condensed polyampholyte dropletis in place, a substantial number of dipole solvent molecules are found inside the droplet.Using methods described above and in previous studies, we determine the phasediagrams for the three polyampholyte sequences sv15, sv20, and sv30, each using five differ-ent solvent models (Fig. 6). The upper critical solution temperature ( T ∗ , cr )—above whichphase separation is impossible—is estimated for each system (Table 1). When comparisonis made across different sequences, LLPS propensity is seen to increase from sv15 to sv20to sv30 for all five solvent models, as quantified by a monotonic increase in T ∗ , cr from leftto right for every row in Table 1, consistent with and generalizing previous findings thatLLPS propensity tend to increase with increasingly negative values of the SCD chargepattern parameter. (a) sv15 (b) sv20 (c) sv30 ,,,, Fig. 4:
Density and dielectric profiles of phase-separated polyampholytes in dipole solvents. Simulationresults for systems each containing n p = 100 polyampholyte chains and n w = 12 ,
000 dipoles are shownfor sequences sv15 (a), sv20 (b), and sv30 (c) at temperatures (top to bottom) T ∗ = 3 .
0, 4 .
0, 5 .
0, and5 . T ∗ = 3 .
0, 4 .
0, 5 .
0, and 6 . T ∗ = 3 .
0, 5 .
0, 6 .
0, and 6 . ρ w , blue curves) and polyampholytes ( ρ , red curves) as a function of the horizontal coordinate z of the simulation box are given by the vertical scale on the right, with ρ w ≈ .
252 outside of thepolyampholyte droplet (peak region of the red curve) in every case. Local effective relative permittivity (cid:15) r ( z ) arising from the dipole molecules (green curves, vertical scale on the left) is calculated by usingEq. 16 for (cid:15) r (cid:107) ( z ). The ρ w ( z ), ρ ( z ), and (cid:15) r ( z ) profiles are computed by averaging over successive bins ofwidth 10 a and plotting the averages at the midpoint of each bin. The effective (cid:15) r for bulk dipole solventis maintained at an essentially constant value of ≈ . T ∗ s for all the systems consideredin this figure by using the q d values in Fig. 3b and q d = ± .
191 for T ∗ = 5 . q d = ± .
618 for T ∗ = 6 . As far as solvent model is concerned, Fig. 6 compares three explicit-solvent models [(i),(ii), and (iii)] and two implicit-solvent models [(iv) and (v)]. The implicit-solvent modelresults for sv20 are simulated here anew, those for sv15 and sv30 are taken from ref. 55.Between these two models, the lower LLPS propensity of model (v) than model (iv)—forall three sequences—follows directly from the definition that the LLPS-driving effectiveelectrostatic interactions among the polyampholyte chains are weakened by a factor of (cid:15) r = 4 . (a)(b) (d)(e)(f)(c) Fig. 5:
Polyampholyte phase separation in explicit dipole solvent. Snapshots of n p = 100 sv20polyampholyte chains and n w = 12 ,
000 dipoles depict a phase-separated condensed polyampholyte dropletat low temperature ( T ∗ = 5 .
5, left column, (a)–(c)) and a lack of phase separation at a higher temperature( T ∗ = 7 .
0, right column, (d)–(f)). Chains and dipoles at the boundaries of the periodic simulation boxesare unwrapped. (a) and (d) show an overview of the entire simulation box, whereas (b) and (e) show across-sectional view of a thin slice parallel to the page. (c) and (f) are enlarged views of the left halvesof (a) and (d), respectively. In (a), (b), (d), and (e), polyampholytes are depicted in red, and dipolesolvent molecules are depicted in cyan. In (c) and (f), the positively and negatively charged beads in thepolyampholytes are in blue and red, whereas the positively and negatively charged beads in the dipoles arein cyan and pink, respectively.
The LLPS propensities of the present explicit-solvent models are all substantially higherthan that of the corresponding implicit-solvent models for the same sequence because,by construction, the non-electrostatic LJ-like part of dipole-polyampholyte interactions ispurely repulsive (Eqs. 6 and 7). This amounts to a hydrophobicity-like, general solvo-phobic driving force against solvation of the polyampholytes in the dilute phase and thusfavoring their condensation. While this feature may be viewed as artifactual, as the LJ-likeinteractions of realistic water models are not purely repulsive, the large offset between thecoexistence curves of implicit-solvent and explicit-solvent models in Fig. 6 is inconsequentialto the physical questions at hand. As argued below, for our purpose, it is only necessary tocompare the LLPS properties of the explicit-solvent models.Among the three explicit-solvent models in Fig. 6, the LLPS propensity of model (i) is al-ways higher than that of model (iii) because they share the same nonpolar solvent (“dipoles”with their charges switched off) but LLPS-driving electrostatic interactions among thepolyampholytes are weakened by a factor of (cid:15) r = 4 . Table 1.
Upper Critical Solution Temperature T ∗ , cr for Polyampholyte LLPS SimulatedUsing Different Solvent Models in Fig. 6. Values of the SCD Sequence Charge PatternParameter for the Sequences sv15, sv20, and sv30 are Provided in Parentheses.solvent sv15 sv20 sv30model (SCD = − .
35) (SCD = − .
37) (SCD = − . . . . . . . . . . . . . . . . (a) sv15 (b) sv20 (c) sv30 (i) (i) (i) (ii)(ii) (ii)(iii) (iii) (iii)(iv) (iv)(iv)(v) (v)(v) Fig. 6:
Phase diagrams of polyampholytes in implicit solvents and explicit dipole and nonpolar solvents.Coexistence curves (i), (ii), and (iii) for sequences sv15 (a), sv20 (b), and sv30 (c) are determined bysimulations of systems each consisting of n p = 100 polyampholyte chains and n w = 12 ,
000 dipoles ( | q d | > q d = 0). The charges q d on the dipoles at different temperatures are adjusted byEq. 15 such that the effective relative permittivity arising from the dipoles are essentially constant, at (cid:15) r ≈ .
0, for all temperatures. Corresponding coexistence curves (iv) and (v), determined by implicit-solvent simulations (no explicit solvent molecules) using n p = 500 polyampholyte chains, are included forcomparison. T ∗ is vacuum reduced temperature, ρ is density of polyampholyte monomer beads in units of1 /a . Lines joining computed data points (circles) are merely a guide for the eye. The five different solventmodels [(i)–(v)] are: (i) Explicit nonpolar solvent, i.e., “dipoles” with q d = 0, (cid:15) r = 1 for polyampholyteelectrostatics (green curves). (ii) Explicit dipole solvent, with different | q d | > T ∗ to maintainan essentially T ∗ -independent effective bulk-solvent (cid:15) r ≈ . (cid:15) r = 1 for all electrostatic interactions (bluecurves). (iii) Explicit nonpolar solvent, same as (i) but (cid:15) r = 4 . (cid:15) r = 1 (red curves). (v) Implicit-solvent, (cid:15) r = 4 . T ∗ value of a coexistence curve for a given system is its upper critical solution temperature T ∗ , cr . The T ∗ , cr values for all the systems considered in this figure are provided in Table 1. (cid:15) r ≈ . (cid:15) r = 4 .
0. At a conceptual level, thiscomparison serves to assess the accuracy, or potential pitfall, of the common practice ofusing the bulk-water (cid:15) r in implicit-solvent simulations of LLPS. Comparisonsbetween the blue (ii) and magenta (iii) curves in Fig. 6 indicate that there is, in general,an appreciable difference in predicted LLPS properties between the two approaches. Theobserved differences are not large. It can be very minor for the sequence with a low − SCDvalue (sv15) but are larger for the sequence with a high − SCD value (sv30). The criticaltemperature T ∗ , cr for the model with explicit dipole solvent molecules and full electrostaticinteractions (blue curve) is higher than the corresponding T ∗ , cr for the model with nonpolarsolvent molecules but reduced electrostatic interactions (magenta curve) for sv20 and sv30.However, the reverse is true for sv15, albeit the difference in this case is very small. Takentogether, these observations suggest that although there is a tendency for implicit-solventsimulations of LLPS using bulk-water (cid:15) r to underestimate LLPS propensity, there is noabsolute rule against it giving a small overestimation of LLPS propensity, as that may welldepend on the structural and energetic details of the system. Analytical Theory of Pure Dipole Solvent Systems.
To assess the generality ofthe above observations from molecular dynamics simulations, we have also developed corre-sponding analytical theories that treat both polyampholyte and dipole degrees of freedom,paying particular attention to excluded volume effects of the polyampholyte chains andthe dipoles as model solvent molecules. Analytical formulations are useful for conceptualunderstanding and hypothesis evaluation. As a theoretical tool, they are complementaryto molecular dynamics simulations. Although analytical theories lack the structural andenergetic details of molecular dynamics models, computation needed for numerical solutionsto analytical formulations such as RPA and field-theoretic simulations (FTS) areoften significantly more efficient than molecular dynamics simulations, thus allowing morephysically informative variations in modeling setup to be explored.As for the explicit-chain, explicit-solvent simulation study above, our analytical devel-opment begins with the baseline case of pure dipoles as a model of polar solvents. Considera system of volume Ω containing n w dipoles subject to a local incompressibility condition0to account for dipole excluded volume, ρ w ( r ) = 1 v , (17)where ρ w ( r ) is the number density at the spatial position r and v is the volume assignedto each dipole. (Note that each dipole is counted once by the ρ w here whereas the beaddensity for the dipoles in the molecular dynamics simulations, denoted by the same symbol ρ w , counts each dipole twice because it contains two beads; hence the ρ w value in thepresent analytical theory formulation is equivalent to ( ρ w /
2) in the molecular dynamicsresults reported above). With only dipoles in the system, v − = n w / Ω. We let r i and µ i denote the position and dipole moment of individual dipole i , and model each dipole as asoft sphere following a Gaussian distribution Γ( r ) of width ¯ a . The local number density ρ w ( r ) and dipole moment density m ( r ) (i.e. the polarization) can then be re-expressed intheir respective smeared forms: ρ w ( r ) = n w (cid:88) i =1 Γ( r − r i ) (18a) m ( r ) = n w (cid:88) i =1 µ i Γ( r − r i ) (18b)with the Gaussian smearing function Γ( r ) = e − r / a / (2 π ¯ a ) / . This smearing procedurewas introduced in ref. 111 and is a convenient way of handling ultraviolet (short-distance)divergences arising from point interactions. A spatially varying m ( r ) gives rise to a boundcharge density ρ c ( r ) = − ∇ · m ( r ). We consider two distinct cases where the dipoles areeither given freely orientable dipole moments of fixed, i.e., permanent, magnitude, | µ i | = µ (“fixed dipoles”), or where every dipole is associated with a harmonic energy contribution κ d µ i / κ d is the corresponding spring constant.The canonical partition function for this system is Z = 1 n w ! n w (cid:89) i =1 (cid:90) d r i (cid:90) d µ i e − U δ (cid:20) ρ w ( r ) − v (cid:21) , (19)where the functional δ -function implements the incompressibility condition in Eq. 17 and U is the Hamiltonian. In the case of fixed dipoles, the integrations over µ i are reduced tosolid angle integrals, and U is the total electrostatic potential energy, U = 12Ω (cid:88) k (cid:54) =0 ˆ ρ c ( k ) 4 πl B k ˆ ρ c ( − k ) = l B (cid:90) d r (cid:90) d r (cid:48) ρ c ( r ) ρ c ( r (cid:48) ) | r − r (cid:48) | (20)1where ˆ ρ c ( k ) = i k · ˆ m ( k ) is the Fourier transform of ρ c ( r ), l B is the vacuum Bjerrum lengthas defined above, and i = −
1. Note that the k = mode can be dropped in the abovesummation because overall charge neutrality of the system entails ˆ ρ c ( k = ) = 0. In the caseof polarizable dipoles, U contains the additional term (cid:80) n w i =1 βκ d µ i / β ≡ /k B T .A formula for the relative permittivity (cid:15) r arising from the dipoles can be obtained byimagining adding a small test charge to the system and looking at the induced electrostaticpotential far away from the test charge. As derived in ref. 112, this consideration yields arelation between (cid:15) r and the charge-charge correlation function:1 (cid:15) r = lim k → (cid:18) − πl B k Ω (cid:104) ˆ ρ c ( k ) ˆ ρ c ( − k ) (cid:105) (cid:19) , (21)where (cid:104) . . . (cid:105) denotes the thermal average over { r i , µ i } .Following a standard procedure using Hubbard-Stratonovich transformations, thepartition function in Eq. 19 can be turned into a statistical field theory where the { r i , µ i } variables are traded in favor of charge- and number density conjugate fields ψ ( r ) and w ( r ),leading to Z = Ω n w n w ! (cid:90) D w D ψ e − H (22)up to an inconsequential overall multiplicative constant, where H [ w, ψ ] = − n w ln Q w + 12Ω (cid:88) k (cid:54) = k πl B ˆ ψ ( k ) ˆ ψ ( − k ) − i v ˆ w ( ) , (23)and (cid:90) D w D ψ = (cid:89) k (cid:54) = (cid:115) k πl B v (cid:90) d ˆ w ( k ) (cid:90) d ˆ ψ ( k ) (24)in k -space representation, with ˆ ψ ( k ) and ˆ w ( k ) being the Fourier transforms, respectively,of ψ ( r ) and w ( r ). The field Hamiltonian H is thus expressed in terms of the Fouriertransformed fields ˆ w and ˆ ψ . In Eq. 23, Q w is a single dipole partition function, whose formfollows from the Hubbard-Stratonovich transformation. It is Q w = 1Ω (cid:90) d r exp (cid:40) − iΩ (cid:88) k ˆ w ( k )ˆΓ( − k )e − i k · r (cid:41) sinh (cid:104) (1 / Ω) (cid:80) k (cid:54) = | k | µ ˆ ψ ( k )ˆΓ( − k )e − i k · r (cid:105) (1 / Ω) (cid:80) k (cid:54) = | k | µ ˆ ψ ( k )ˆΓ( − k )e − i k · r (25)2for fixed dipoles, and Q w = 1Ω (cid:90) d r exp (cid:40) − iΩ (cid:88) k ˆ w ( k )ˆΓ( − k )e − i k · r (cid:41) exp (cid:104)(cid:80) k (cid:54) = k ˆ ψ ( k )ˆΓ( − k )e − i k · r (cid:105) βκ d Ω (26)for polarizable dipoles, where the Fourier-transformed smearing function ˆΓ( k ) = e − ¯ a k / .Recognizing that the charge-charge correlation function in Eq. 21 may be rewritten in afield-theory representation involving the correlation function of the ψ field, (cid:104) ˆ ρ c ( k ) ˆ ρ c ( − k ) (cid:105) = k πl B (cid:34) − k πl B (cid:104) ˆ ψ ( k ) ˆ ψ ( − k ) (cid:105) Ω (cid:35) , (27)and substituting this expression into Eq. 21 leads to the field-theory expression for therelative permittivity, 1 (cid:15) r = lim k → k πl B Ω (cid:104) ˆ ψ ( − k ) ˆ ψ ( k ) (cid:105) , (28)where now (cid:104) . . . (cid:105) denotes an average over field configurations weighted by e − H .We next proceed to the calculation of this ψ correlation function in order to compute (cid:15) r . As a first approximation—which corresponds to RPA—we compute it analytically ac-counting only for Gaussian fluctuations in the fields. Introducing | Ψ( k ) (cid:105) = [ ˆ w ( k ) , ˆ ψ ( k )] T ,where the superscript “T” denotes (vector) transposition, the field Hamiltonian in thisapproximation is H ≈ (cid:88) k (cid:54) = (cid:104) Ψ( − k ) | ˆ∆( k ) | Ψ( k ) (cid:105) , (29)which contains only terms quadratic in the fields, and whereˆ∆( k ) → ˆ∆ k = (cid:32) ˆΓ k /v
00 ( k / πl B ) (cid:16) χ D ˆΓ k (cid:17)(cid:33) (30)with k = | k | , ˆΓ k standing for e − ¯ a k / , and χ D = πl B µ v (fixed dipoles)4 πKv (polarizable dipoles) (31)where K ≡ l B /βκ d is a temperature-independent factor corresponding to the molecularpolarizability in our reduced units. When Eq. 28 is computed using the RPA HamiltonianEq. 29, the ψ correlation is simply the ψψ element of ˆ∆ − , which is equal to the reciprocal3of the ψψ component of ˆ∆ (bottom-right element) in Eq. 30. It follows that (cid:15) r = 1 + χ D (32)in RPA, yielding the standard formula in classical electrodynamics for uncorrelated dipoles.Eq. 32 with the χ D expression for fixed dipoles in Eq. 31 is equivalent to Eq. 14 for analyzingmolecular dynamics simulation data on pure dipole systems in Fig. 2.Going beyond RPA, we can calculate the ψ correlation function more accurately by usingFTS on a lattice. In FTS, the continuum fields w ( r ) and ψ ( r ) are approximated by discretefield variables defined on the sites of a cubic lattice with periodic boundary conditions.Due to the complex nature of H , field averages (such as the ψ correlation function) aretypically computed in FTS using the Complex-Langevin (CL) method, —a techniqueinspired by stocahstic quantization of field theories —where a fictitious time coordinateis introduced such that w ( r ) → w ( r , t ) and ψ ( r ) → ψ ( r , t ) and the fields are now complex.The fields evolve in CL time according to the Langevin equation ∂ϕ ( r , t ) ∂t = − δ H δϕ ( r , t ) + η ϕ ( r , t ) , ϕ = w, ψ , (33)where η ϕ is a real-valued Gaussian noise with zero mean and (cid:104) η ϕ ( r , t ) η ϕ (cid:48) ( r (cid:48) , t (cid:48) ) (cid:105) = 2 δ ϕ,ϕ (cid:48) δ ( r − r (cid:48) ) δ ( t − t (cid:48) ). In the CL method, field averages can be shown to correspond to asymptotic CLtime averages. In this work, we use the first order semi-implicit method of ref. 121 to solveEq. 33 numerically. l B (cid:15) r − − T ∗ fi x e d d i p o l e s , µ = . fi x e d d i p o l e s , µ = . polarizable dipoles Fig. 7:
Relative permittivity computed by FTS for incompressible pure dipoles. Each error barindicates the standard deviation from three independent FTS runs. Dashed lines correspond to RPAEq. 32. The reduced temperature, T ∗ (top horizontal scale), used in the analytical theories presentedin this work is defined by T ∗ = l/l B , where l = √ a is the reference polymer bond length of thepolyampholyte chains to be discussed below. The magnitude of the dipole moment, µ , is in unitsof el . l is equivalent to the reference bond length a used in Eq. 12 for T ∗ in the molecular dynamics study. (cid:15) r computed by FTS for several representative parameter values. The FTS ψ correlation functions were computed on a 24 lattice with a lattice spacing ∆ x = ¯ a = 1 / √ t = 10 − . After an initial equilibration period of 2 × time steps, thefields were sampled every 1,000 steps for ∼ × steps. Fig. 7 shows the mean and stan-dard deviation among three independent simulations that were performed for each valueof l B . All simulations used v = 0 .
1. The results for fixed dipoles were obtained for µ = 0 . µ = 0 .
35. The polarizable dipoles were simulated using K = 0 . (cid:15) r = 20 according to Eqs. 31 and 32. With the relative permittivity of the polymer materialset to unity in the following analytical theory development, the temperature-independent (cid:15) r value of 20 corresponds approximately to the ratio of the relative permittivities of ≈ ∼ Analytical Theory for Mixture of Dipole and Nonpolar Solvents.
Before intro-ducing polymers to the pure dipole systems described above, it is instructive to first considerthe simpler case of an incompressible mixture of fixed dipoles and nonpolar (and nonpo-larizable) particles as a model of nonpolar solvent molecules. For this purpose, we add n n Gaussian-smeared particles with positions r n ,i and modify the incompressibility conditionin Eq. 17 to ρ w ( r ) + ρ n ( r ) = 1 v , (34)where ρ n ( r ) = (cid:80) n n i =1 Γ( r − r n ,i ). The canonical partition function for this system is Z = 1 n w ! n n ! n n (cid:89) i =1 (cid:90) d r n ,i n w (cid:89) i =1 (cid:90) d r i (cid:90) d µ i e − U δ (cid:20) ρ w ( r ) + ρ n ( r ) − v (cid:21) , (35)where U is unchanged. Accordingly, the field theory version of the partition function isgiven by Z = Ω n w + n n n w ! n n ! (cid:90) D w (cid:90) D ψ e − H , (36)where now H [ w, ψ ] = − n w ln Q w − n n ln Q n + 12Ω (cid:88) k (cid:54) = k πl B ˆ ψ ( k ) ˆ ψ ( − k ) − i v ˆ w ( ) (37)includes a term with the partition function Q n for a single nonpolar particle, viz., Q n = 1Ω (cid:90) d r exp (cid:40) − iΩ (cid:88) k ˆ w ( k )ˆΓ( − k )e − i k · r (cid:41) . (38)5With this setup, once again we can apply Eq. 28 to compute the bulk relative permittivityof mixtures of the two solvents, both in RPA and FTS. In RPA, we have that (cid:15) r = 1 + χ D φ w . (39)where φ w = vn w / Ω is the dipole volume fraction. The corresponding FTS calculation fol-lows the same recipe as under the previous subheading, but now with the field Hamiltonianin Eq. 37. Comparisons between (cid:15) r computed in RPA and FTS are shown in Figs. 8 and 9for µ = 0 . µ = 0 .
35, respectively, for a wide range of l B values. Common to all runsare the parameter values v = 0 . a = 1 / √
6. FTS follows RPA quite well except athigh l B and intermediate φ w where the FTS (cid:15) r can be substantially lower than the RPA (cid:15) r (see, e.g., Fig. 9c). (cid:15) r (a) l B = 10 (b) l B = 100 − − φ w (cid:15) r (c) l B = 200 − − φ w (d) l B = 300 Fig. 8:
Relative permittivity (cid:15) r computed in FTS and RPA for an incompressible mixture of fixed dipolesand nonpolar particles with µ = 0 . l B (in units of l ). As in Fig. 7, error bars arestandard deviations of independent FTS runs, and dashed curves are RPA predictions. − − φ w (cid:15) r (a) l B = 10 − − φ w (b) l B = 70 − − φ w (c) l B = 150 Fig. 9:
Same as Fig. 8 but for µ = 0 . (cid:15) r may be understood by envisioning the system entering into an inho-mogeneous state where a dipole-rich and a dipole-depleted phase coexist. Obviously, byconstruction, this phase separation is entirely driven by the electrostatic dipole-dipole in-teractions; and the system may be viewed as a special case of a binary mixture of twotypes of dipoles distinguished only by their dipole moments, which has been considered inref. 92. In this reference, the virial expansion of such a system was studied using RPA, andit was shown that there is an effective Flory χ parameter, proportional to the square of thedifference in dipole moment, that can induce phase separation. By directly calculating Z inRPA, we can subsequently compute the binodal curve using the free energy f = − ln Z / Ω.In RPA, the field Hamiltonian still follows Eq. 29, but now withˆ∆ k = (cid:32) ˆΓ k /v
00 ( k / πl B ) (cid:16) χ D ˆΓ k φ w (cid:17)(cid:33) , (40)which leads to the free energy f = φ w ln φ w + (1 − φ w ) ln(1 − φ w ) + v π (cid:90) ∞ d k k ln (cid:104) χ D φ w ˆΓ (cid:105) , (41)up to an additive constant. Fig. 10 shows the phase diagram computed for the two values ofthe dipole moment in Figs. 8 and 9 using this free energy. Dotted horizontal lines indicatethe l B values considered in FTS. Indeed, we see that the parameter values for which FTSand RPA do not quite agree with regard to (cid:15) r ( l B = 300 for µ = 0 . l B = 150for µ = 0 .
35 in Fig. 9) coincide well with parameter values in Fig. 10 for which RPApredicts phase separation with significant differences in dilute- and condensed-phase φ w s. − − − − φ w l B µ = 0 . µ = 0 . Fig. 10:
Phase diagrams computed in RPA for incompressible mixtures of fixed dipoles and nonpolarparticles. The binodal (coexistence) curves are shown for two values of the dipole moment, µ = 0 . µ = 0 . v = 0 .
1, and ¯ a = 1 / √ l , as noted for Fig. 7). The dashed horizontal lines indicate l B values considered in FTS. Analytical Theory for Polyampholytes in Dipole Solvent.
We now consider sys-tems consisting of n w dipoles and, instead of nonpolar particles, n p polyampholytes. Fol-lowing the notation in ref. 51, each polyampholyte (labeled by α ) is made up of N monomerbeads (as in the above molecular dynamics simulations) at positions R α,τ , τ = 1 , , . . . , N ,and carries charge sequence | σ (cid:105) = [ σ , σ , ...σ N ] T . The incompressibility constraint for thissystem is given by ρ p ( r ) + ρ w ( r ) = 1 v (42)wherein the polymer matter (mass) density ρ p ( r ) = (cid:80) n p α =1 (cid:80) Nτ =1 Γ( r − R α,τ ). The Fourier-transformed density of polymer charges,ˆ c ( k ) = ˆΓ( k ) n p (cid:88) α =1 N (cid:88) τ =1 σ τ e i k · R α,τ (43)now contribute to the Fourier-transformed total charge density such thatˆ ρ c ( k ) = ˆ c ( k ) + i k · ˆ m ( k ) . (44)The partition function of the polyampholytes-plus-dipoles system is given by Z = 1 n p ! n w ! n p (cid:89) α =1 N (cid:89) τ =1 (cid:90) d R α,τ n w (cid:89) i =1 (cid:90) d r i (cid:90) d µ i e − T [ R ] − U [ R , r ] δ (cid:20) ρ w ( r ) + ρ p ( r ) − v (cid:21) , (45)where T [ R ] = 32 l n p (cid:88) α =1 N − (cid:88) τ =1 ( R α,τ +1 − R α,τ ) (46)is the Hamiltonian term for Gaussian chain connectivity, where l is the reference bondlength that we have used as a length scale in Fig. 7 to define T ∗ for the present analyticaltheories, and U for electrostatic interactions takes the same general form as that given byEq. 20. The field-theory expression for the partition in Eq. 45 is Z = Ω n p + n w n p ! n w ! (cid:90) D w D ψ e − H (47)where H [ ψ, w ] = − n p ln Q p − n w ln Q w + 12Ω (cid:88) k (cid:54) = k πl B ˆ ψ ( k ) ˆ ψ ( − k ) − i v ˆ w ( ) (48)and, via Hubbard-Stratonovich transformation, the single-polymer partition function Q p is8given by Q p = 1Ω N (cid:89) τ =1 (cid:90) d R τ exp (cid:40) − l N − (cid:88) τ =1 ( R τ +1 − R τ ) − iΩ (cid:88) k N (cid:88) τ =1 (cid:104) ˆ w ( k ) + σ τ ˆ ψ ( k ) (cid:105) ˆΓ( − k )e − i k · R τ (cid:41) . (49)As in the above analysis of pure dipoles and dipoles-plus-nonpolar-particles systems, wemay now expand the field Hamiltonian to second order in the fields to arrive at an RPAtheory. Again, this leads to the general RPA form in Eq. 29, now with a different expressionfor the kernel ˆ∆ k for the system described by Eq. 45:ˆ∆ k = (cid:32) ( φg mm k + φ w )ˆΓ k /v φg mc k ˆΓ k /vφg mc k ˆΓ k /v ( φg cc k + D φ w k )ˆΓ k /v + ( k / πl B ) (cid:33) , (50)where D = µ / D = 1 / ( βκ d ) for induced (polarizable) dipoles.Here φ = vN n p / Ω is the polymer volume fraction (thus φ w = 1 − φ is the dipole volumefraction). The factors g mm k , g mc k and g cc k stem from the single-polymer partition function inEq. 49, g mm k = 1 N N (cid:88) ν,τ =1 exp( − l k | ν − τ | / , (51a) g mc k = 1 N N (cid:88) ν,τ =1 σ ν exp( − l k | ν − τ | / , (51b) g cc k = 1 N N (cid:88) ν,τ =1 σ ν σ τ exp( − l k | ν − τ | / . (51c)Applying the RPA form of the Hamiltonian H in Eq. 29 to the expression for Z in Eq. 22and performing the (cid:82) D w D ψ integrals by including the prefactors of (cid:82) d ˆ w ( k ) (cid:82) d ˆ ψ ( k ) inEq. 24 lead to the following expression for the free energy f = − ln Z / Ω in RPA: f = φN ln φ + (1 − φ ) ln(1 − φ ) + v π (cid:90) ∞ d k k ln (cid:20) πl B vk det ˆ∆ k (cid:21) . (52)where ˆ∆ k is here given by Eq. 50. We rewrite the argument of the logarithm in the aboveintegrand in the form of (4 πl B v/k ) det ˆ∆ k (4 πl B v/k ) det ˆ∆ k (cid:12)(cid:12)(cid:12) φ =0 = 1 + A k φ + B k φ (53)9to arrive at a form of f in which irrelevant constant terms are subtracted, with A k = g mm k − πl B vk (1 + χ D ˆΓ k ) ˆΓ k ( g cc k − k D ) , (54a) B k = 4 πl B vk (1 + χ D ˆΓ k ) ˆΓ k (cid:2) ( g mm k − g cc k − k D ) − ( g mc k ) (cid:3) . (54b)At this juncture, we are in a position to use Eq. 52 to compute phase diagrams in RPA.Fig. 11 shows such phase diagrams in the ( φ, T ∗ ) plane for polyampholyte chains with chargesequence sv15 immersed in fixed dipoles (Fig. 11a) and polarizable dipoles (Fig. 11c). Notsurprisingly, compared to the case with nonpolar solvent (solid curves), LLPS propensityis lowered in the presence of polar solvents (dashed, dotted, and dashed-dotted curvesin Figs. 11a and c) because the LLPS-driving effective electrostatic interactions amongthe polyampholytes are weakened due to screening, and this LLPS-suppressing effectincreases with increasing polarity of the solvent ( T ∗ , cr lower for µ = 0 . µ = 0 . φ, T ∗ ) plane, where T ∗ = T ∗ (cid:15) r ( φ = 0) with (cid:15) r ( φ = 0) being therelative permittivity of the bulk-solvent. These plots offer a comparison between the phasebehavior of polyampholytes interacting via a homogeneous background dielectric constant(solid curves) and the phase behavior of polyampholytes interacting with an incompressibledipole solvent (thus allowing for dipole-concentration-dependent relative permittivity) withthe same bulk-solvent relative permittivity. Thus this comparison corresponds to thatbetween molecular dynamics models (iii) and (ii) in Fig. 6a. Interestingly, whereas models(iii) and (ii) in Fig. 6a exhibit essentially the same LLPS propensity, RPA predictions inFigs. 11b and d indicate that interacting with an incompressible dipole solvent enhancespolyampholyte LLPS propensity relative to that predicted for the corresponding situationwith a homogeneous background dielectric constant. This enhancing effect is quite smallfor the fixed dipoles in Fig. 11b—the effect nonetheless increases with µ —but the effect ismore appreciable for the case of polarizable dipoles in Fig. 11d. As discussed previously, because of the relatively strong attractive LJ background in the molecular dynamicsmodel—but this feature is not considered in the present theoretical formulation, theupward concavity on the condensed sides of the RPA coexistence curves in Fig. 11 is notobserved for the molecular dynamics coexistence curves in Fig. 6.0 T ∗ (a) µ = 0 . µ = 0 . fixed dipoles (c) K = 0 . polarizable dipoles . . . φ T ∗ (b) 0 . . . φ (d) Fig. 11:
RPA phase diagrams of incompressible mixtures of polyampholytes and dipoles. Results are forsequence sv15 with v = 0 . a = 1 / √ l ). φ is polymer volume fraction, T ∗ = l/l B is vacuum reduced temperature as in Fig. 7, and T ∗ = T ∗ (cid:15) r ( φ = 0) is reduced temperature that takesinto account bulk-solvent relative permittivity (which is temperature dependent for fixed dipoles). Severalsolvent models are compared. Solid curves (all panels) are for nonpolar solvent ( µ = 0); dashed and dottedcurves [(a) and (b)] are for fixed-dipole solvents with µ = 0 . µ = 0 .
1, respectively; and dashed-dottedcurves [(c) and (d)] are for polarizable-dipole solvent with K = 0 . (cid:15) r ( φ = 0) = 20). Note that the l B values considered in this figureare much lower than that required for the dipoles to phase separate on their own (Figs. 8–10). Comparison of Theory-Predicted Polyampholyte Phase Behaviors Under Dif-ferent Relative Permittivity Models.
To clarify the difference between various per-mittivity models that have been employed by our group to study polyampholyte LLPSusing analytical theories, we take the mean-field limit for mass-density fluctuation,i.e. w →
0, of the above formulation to focus on features of the RPA electrostatic freeenergy. In such a limit, only the bottom-right component of the ˆ∆ array in Eq. 50 remains,viz. ˆ∆ k −→ w → v (cid:2) φg cc k + (1 − φ ) D k (cid:3) ˆΓ k + k πl B . (55)1The RPA electrostatic free energy is then given by f el = v π (cid:90) ∞ d k k ln (cid:20) χ D ˆΓ k (1 − φ ) + 4 πl B vk g cc k ˆΓ k φ (cid:21) , (56)which may be rewritten as f el = v π (cid:90) ∞ d k k ln (cid:104) χ D ˆΓ k (1 − φ ) (cid:105) + v π (cid:90) ∞ d k k ln (cid:34) πl B g cc k ˆΓ k φvk [1 + χ D ˆΓ k (1 − φ )] (cid:35) ≡ f ddel + f (cid:48) el (57)in which the first term f ddel accounts for dipole-dipole interaction between solvent molecules,and the second term includes charge-charge (polymer-polymer) and charge-dipole (polymer-solvent) interactions. It should be noted that in the present system, the pure material sus-ceptibilities are χ pD = 0 for polymers (polyampholytes) and χ wD = χ D for solvent molecules(dipoles). A concentration-dependent permittivity which for analytical tractability is as-sumed to be dependent linearly on volume fractions can then be constructed as (cid:15) r ( φ ) = (cid:15) p φ + (cid:15) w (1 − φ )=(1 + χ pD ) φ + (1 + χ wD )(1 − φ )= φ + (1 + χ D )(1 − φ )=1 + χ D (1 − φ ) . (58)Thus, in the case of vanishing Gaussian smearing, i.e. ¯ a → k →
1, the second termin Eq. 57 can be cast in the following form:lim ˆΓ k → f (cid:48) el = f empel = v π (cid:90) ∞ d k k ln (cid:20) πl B v(cid:15) r ( φ ) k g cc k φ (cid:21) , (59)where the superscript “emp” stands for “empirical permittivity”. It is noteworthy that f empel derived above is equivalent to the RPA electrostatic energy of the recent “with ‘self-energy’” model used in Fig. 7B of Das et al. More specifically, with v = l , f empel is seento be identical to the expression in Eq. S60 in the SI Appendix of ref. 34, which in turn isbased on the partition function Z el in Eq. S51 of the same reference.We now proceed to clarify how dipole-dipole interactions promote LLPS whereas dipole-polymer interactions mostly suppress LLPS in the present analytical models. First, a small- φ Taylor expansion of the logarithmic integrand of f ddel in Eq. 57 yieldsln (cid:104) χ D ˆΓ k (1 − φ ) (cid:105) = ln (cid:104) χ D ˆΓ k (cid:105) − χ D ˆΓ k χ D ˆΓ k φ − (cid:32) χ D ˆΓ k χ D ˆΓ k (cid:33) φ + O ( φ ) , (60)2wherein φ is seen to associate with a negative factor, indicating that dipole-dipole inter-actions contribute an effectively attractive polymer-polymer interaction. In other words,dipole-dipole interactions contribute a positive Flory-Huggins χ parameter or a negativesecond virial coefficient and therefore promote polymer phase separation. It should benoted as well that an equivalent conclusion can be reached by considering a Taylor expan-sion with respect to the mean-field concentration of the solvent molecules, (1 − φ ), leadingto the conclusion that the effective interactions among the dipoles are also attractive.Second, we note that the f (cid:48) el term in Eq. 57 for the “empirical permittivity” model maybe dissected further by rewriting it in the following form: f (cid:48) el = v π (cid:90) ∞ d k k ln (cid:20) πl B vk g cc k ˆΓ k φ (cid:21) + v π (cid:90) ∞ d k k ln − πl B vk (cid:16) χ D χ D (cid:17) g cc k ˆΓ k φ πl B vk g cc k ˆΓ k φ . (61)Here the first term is the RPA theory with constant vacuum permittivity, i.e. purely charge-charge (polymer-polymer) interactions, and the second term accounts for charge-dipole(polymer-solvent) interactions, and we have defined χ D ≡ χ D Γ k (1 − φ ). A Taylor expansionof the second term with respect to φ givesln − πl B vk (cid:16) χ D χ D (cid:17) g cc k ˆΓ k φ πl B vk g cc k ˆΓ k φ = − A k B k B k + 1 φ + A k B k ( A k B k + 2 A k + 2)2( B k + 1) φ + O ( φ ) , (62)where A k ≡ πl B vk g cc k ˆΓ k , B k = χ D ˆΓ k . (63)Given A k , B k > ∀ k , the polymer-dipole interactions entail a positive second virial coefficientassociated with the φ term in Eq. 62 that equals to v π (cid:90) ∞ d k k A k B k ( A k B k + 2 A k + 2)2( B k + 1) > , (64)which acts as an effective excluded-volume-like repulsion between polymers and thus sup-press phase separation.The comparison above is between two models with the same polymer relative permittivity( (cid:15) p = 1) but different solvent relative permittivities ( (cid:15) w = 1 v.s. (cid:15) w = 1 + χ D , see Eq. 58).Another comparison, of more direct interest to the question at hand, is between an implicit-solvent model in which the polymers interact in a dielectric environment with a constant (cid:15) w = 1 + χ D (the relative permittivity of the polymer material itself remains the same at (cid:15) p = 1) and the full theory with the polymers immersed in dipole solvent wherein the bulk-dipole relative permittivity is (cid:15) w = 1 + χ D . This comparison may be analyzed by using3a modified Bjerrum length, l rB = l B /(cid:15) w , as characteristic temperature factor, to regroupEq. 61 for f (cid:48) el into a sum of two contributions, f (cid:48) el = v π (cid:90) ∞ d k k ln (cid:20) πl rB vk g cc k ˆΓ k φ (cid:21) + v π (cid:90) ∞ d k k ln πl rB vk g cc k ˆΓ k (cid:16) χ D χ D − (cid:17) φ πl rB vk g cc k ˆΓ k φ , (65)such that the first term is the electrostatic free energy for polymers interacting with aconstant implicit-solvent (cid:15) w = 1 + χ D . The integrand of the second term in Eq. 65, whichis a polymer-solvent interaction term, can now be expanded with respect to φ :ln πl rB vk g cc k ˆΓ k (cid:16) χ D χ D − (cid:17) φ πl rB vk g cc k ˆΓ k φ = − A r k χ D (ˆΓ k − χ D ˆΓ k + 1 φ + A r k χ D (cid:110) A r k (ˆΓ k − (cid:104) χ D (ˆΓ k + 1) + 2 (cid:105) + 2(1 + χ D )ˆΓ k (cid:111) (cid:16) χ D ˆΓ k + 1 (cid:17) φ + O ( φ ) , (66)where A r k = 4 πl rB / ( vk ) g cc k ˆΓ k . Now, the form of the the second virial coefficient, denoted v pw , r2 , that associates with the φ term, v pw , r2 = v π (cid:90) ∞ d k k A r k χ D (cid:110) A r k (ˆΓ k − (cid:104) χ D (ˆΓ k + 1) + 2 (cid:105) + 2(1 + χ D )ˆΓ k (cid:111) (cid:16) χ D ˆΓ k + 1 (cid:17) , (67)becomes rather intricate. To analyze the sign of v pw , r2 , we substitute the definition of A rk into the numerator of the integrand in the above expression and rewrite it as4 πl rB χ D v g cc k k ˆΓ k (cid:20) πl rB vk g cc k (cid:16) ˆΓ k − (cid:17) (cid:16) χ D ˆΓ k + χ D + 2 (cid:17) + 2(1 + χ D ) (cid:21) . (68)For overall neutral polyampholytes, g cc k (cid:28) k → g cc k → k → ∞ , and g cc k does not increase monotonically with k but has a peak at some intermediate k for mostcases (except, e.g., for the strictly alternating sequence sv1, which monotonically increasesfrom 0 to 1 when k increases). For k (cid:28) / ¯ a , ˆΓ k →
1, the general expression in Eq. 68becomes 8 πl rB χ D (1 + χ D ) v g cc k k (cid:20) − πl rB a v g cc k + 1 (cid:21) > , (69)4and for k (cid:29) / ¯ a , ˆΓ k →
0, the same general expression becomes8 πl rB χ D (1 + χ D ) v g cc k k ˆΓ k > . (70)Thus, Eq. 68 contributes a repulsive effect to the integral for v pw , r2 in Eq. 67 in both thesmall- k and large- k regimes. Moreover, the ˆΓ k factor in the same expression tends to reducelarge- k contributions to the v pw , r2 integral. Therefore, it is expected that as long as the peakof g cc k is at k (cid:46) /a , effective polymer-solvent interactions should be repulsive. However,in view of the intricate form of v pw , r2 in Eq. 67, the possibility of effective polymer-solventinteractions being attractive in certain extreme scenarios cannot be precluded.Fig. 12a compares the theory-predicted phase behaviors of sequences sv15, sv20, and sv30under three different relative permittivity models for the solvent. All three sequences studiedexhibit the same ranking ordering of LLPS propensities the three models (as quanitified bytheir critical temperatures):explicit dipole solvent > constant permittivity > empirical permittivity . While recognizing that two of the solvent models investigated in Fig. 12a use the theoreticalformulation with the matter conjugate field w (Eqs. 52–54), it is noteworthy that the rankordering of LLPS propensity in Fig. 12a is consistent with above w → > empiricalpermittivity. Second, after subtracting the free energy for the constant permittivity modelfrom the free energy f (cid:48) el for the empirical permittivity model, a LLPS-disfavoring repulsiveterm probably remains (Eqs. 65, 67–70 and discussion that follows), i.e., it is likely thatconstant permittivity > empirical permittivity. These analytical trends thus provide a par-tial rationalization of the rank ordering of LLPS propensity in Fig. 12a, although the above w → While the procedure isvalid for polyampholytes interacting in a medium of uniform relative permittivity because inthat case the self-energy term does not affect phase behaviors, the procedure is problematicfor systems with polymer-concentration-dependent relative permittivity because in thosesituations the self-energy term may capture part of the polymer-solvent interactions andtherefore can impact phase equilibrium, as we have pointed out recently. The presentanalysis allows for further evaluation of such a self-energy term. Consider, in the theoretical5framework developed above, the expression f self ≡ (cid:90) ∞ d kk π πl B (cid:15) r ( φ ) k N ˆΓ k φ N (cid:88) τ =1 σ τ , (71)which can be seen to be equivalent to the self-energy G (˜ k ) quantity in Eq. 69b of ref. 28when there is no salt or counterion ( φ s = φ c = 0), no short-range cutoff of Coulombinteractions (˜ k [1 + ˜ k ] → k ), and η = 1, by noting that k = ˜ k/l ( l now corresponds to b = a in ref. 28), l B = l/T ∗ , ˆΓ k → a →
0, and that σ i = | σ i | for polyampholytes with σ i = ±
1. An expansion of the factor [ (cid:15) r ( φ )] − = [1 + χ D (1 − φ )] − (Eq. 58) in the integrandfor f self in Eq. 71 yields f self = (cid:90) ∞ d kk π πl B k N ˆΓ k φ N (cid:88) τ =1 σ τ (cid:2) − χ D (1 − φ ) + O (cid:0) χ (cid:1)(cid:3) . (72)This form conveys two messages. First, while a part of f self may be identified with a partof f (cid:48) el , e.g., the term linear in φ in the expansion of f self is equal to a part of the first term(linear in φ ) in the expansion of the logarithmic integrand of the first integral in Eq. 61(the µ = ν part of the summation for g cc k in Eq. 51c), and therefore may be viewed ascapturing part of polymer-solvent interactions, f self by itself does not coincide with thefree energy arising from polymer-solvent interactions. Second, the φ term in the f self expanion in Eq. 72 is positive, indicating that f self entails an effective polymer-polymerrepulsion. It follows that subtracting f self should increase LLPS propensity. This effectis illustrated by the coexistence curve in Fig. 12b computed using a free energy with f self subtracted (dotted red curve). The dramatic increase in LLPS propensity of thismodel compared with corresponding models without the self-energy subtraction (orangecurves) is in line with some of our previous results on models with self-energy subtraction(Fig. 12 of ref. 28). However, for other heteropolymer models with less charge densities andaugmented Flory-Huggins interactions, subtracting an electrostatic self-energy term canlead to more modest—though still substantial—increases in LLPS propensity (Fig. 7C ofref. 34). Most importantly, the present analysis shows that a self-energy term similar to theone in Eq. 71 does not account neatly for polymer-solvent interactions, and in any event,subtraction of such a term from the RPA electrostatic free energy is physically unwarranted.6 (a) (b)sv15sv20 sv30 Fig. 12:
Comparing RPA LLPS theories with different solvent permittivity models. The vertical variable T ∗ (cid:15) r = T ∗ , where (cid:15) r is bulk-solvent relative permittivity. (a) Coexistence curves of sequences sv15, sv20,and sv30 are shown in different colors. Dashed-dotted curves are for the full explicit dipole solvent model,solid curves are for the implicit-solvent constant permittivity model (with nonpolar solvent). Both set ofcurves are computed with the mass density conjugate field w in accordance with the formulation given byEqs. 52–54 (dashed-dotted and solid curves for sv15 are the same as those in Fig. 11d). Dashed curvesare for the empirical permittivity model, computed in the w → v = 0 . l ). (b) Coexistencecurves for sv20 with the same models in (a) depicted in the same line styles in orange but all computedusing the w → v = 1 . l B → l rB = l B /(cid:15) w for dashed-dotted curve, first termin Eq. 65 for solid curve, and second term in Eq. 57, as in (a), for dashed curve). The red dotted curve fora “no self energy” model is the coexistence curve computed using an electrostatic free energy equals tothat of the empirical permittivity model (second term in Eq. 57) minus the “self-energy” term in Eq. 71. DISCUSSION
FTS of Polyampholytes with Dipole Solvent.
In addition to RPA, we have ex-tended our FTS effort beyond the study of pure dipoles and dipole and nonpolar solventmixtures (Figs. 7–9) to the investigation of polyampholyte LLPS in the presence of dipoles.In this endeavor, we find that the strict incompressibility constraint considered above forRPA (Eq. 42) is challenging to realize numerically in the fully fluctuating FTS simulationsof systems comprised of polyampholytes and explicit polar solvents. Therefore, we optinstead for a compressible model in which the Dirac δ [ ρ p ( r ) + ρ w ( r ) − /v ] for incompress-ibility is replaced by an energy penalty term in the interaction Hamiltonian in the form of (cid:82) d r ( ρ p ( r ) + ρ w ( r ) − /v ) / γ . The resulting system is incompressible in the γ → γ . The corresponding field-picture Hamiltonian is given by H [ ψ, w ] = 18 πl B (cid:90) d r ( ∇ ψ ( r )) + γ (cid:90) d r w ( r ) − i v (cid:90) d r w ( r ) − n w ln Q w − n p ln Q p , (73)where the partition functions Q w and Q p for a single dipole and a single polymer chain,respectively, are given by Eq. 25 and Eq. 49. We focus here on a solution comprised ofa single species of polyampholytes and a single species of dipoles as a model of a polarsolvent. The equilibrium averages of thermodynamic variables of the system are obtainedby averaging over sufficiently large number of statistically independent configurations ofthe fields ψ and w . Following the CL prescription in Eq. 33, the field configurations aregenerated by numerically solving the following two coupled differential equations: ∂w ( r ) ∂t = − (cid:20) i (cid:18) ˜ ρ p ( r ) + ˜ ρ w ( r ) − v (cid:19) + γw ( r ) (cid:21) + η w ,∂ψ ( r ) ∂t = − (cid:20) i (˜ c ( r ) + ˜ c w ( r )) − πl B ∇ ψ ( r ) (cid:21) + η ψ . (74)In Eq. 74, t is a fictitious time, η w and η ψ are real-valued random numbers drawn from aNormal distribution with zero mean and variance 2 δ ( t − t (cid:48) ) δ ( r − r (cid:48) ), and the field operatorsfor polymer bead density ( ˜ ρ p ), polymer charge density (˜ c ), solvent bead density ( ˜ ρ w ) andsolvent charge density (˜ c w ) are given by˜ ρ p ( r ) = i n p Q p Ω δ Q p δw , ˜ c ( r ) = i n p Q p Ω δ Q p δψ , ˜ ρ w ( r ) = i n w Q w Ω δ Q w δw , ˜ c w ( r ) = i n w Q w Ω δ Q w δψ . (75)The individual solvent beads are overall charge neutral, but polar in nature with eachhaving a fixed magnitude of dipole moment equals to µ . ˜ c w can be re-expressed by defininga solvent polarization density operator such that ˜ c w = − ∇ · (cid:101) m ( r ) where (cid:101) m ( r ) = i n w Q w Ω e − i ˘ w ( r ) ∇ ˘ ψ ( r ) (cid:104) ∇ ˘ ψ ( r ) (cid:105) cos (cid:16) µ | ∇ ˘ ψ ( r ) | (cid:17) − sin (cid:16) µ | ∇ ˘ ψ ( r ) | (cid:17) µ | ∇ ˘ ψ ( r ) | , (76)and the smeared fields in this expression, each marked by a diacritic breve, are defined byspatial convolution as ˘ ϕ = Γ (cid:63) ϕ , ϕ = w, ψ . We integrate Eq. 74 numerically by usingthe first-order semi-implicit time integration method of ref. 121 over a cubic box of volumeΩ = (32¯ a ) . To achieve good numerical precision, we discretize the box into a 32 × × v = 0 .
1, the bulk polymer bead densityis set at 0 . µ = 0 . γ = 10. Although not rigorously identical to matter or charge8densities, the CL time snapshots of density operators are an intuitive and instructive wayto visualize phase separation in FTS systems. Fig. 13 shows representative snapshotsof the field operators of different observables for two different temperatures to illustratethat the present polyampholytes plus dipoles system can indeed undergo phase separation.As temperature is lowered from a relatively high temperature ( T ∗ = 50) to a relatively lowtemperature ( T ∗ = 0 . (a) (b) (c)(d) (e) T ∗ = 0 . T ∗ = 50 (f) Fig. 13:
Illustrative FTS snapshots of the density operators at low and high temperatures. All resultsshown are for polyampholyte sequence sv20, dipole µ = 0 .
1, excluded-volume parameter v = 0 .
1, and γ = 10. (a) Temperature T ∗ = 0 .
2. Real non-zero parts of ˜ ρ p ( r ) and ˜ ρ w ( r ). A droplet-like structure ofsv20 molecules (shown in red) is buried inside the polar solvent (shown in cyan). (b) A cut-plane view of(a) of a thin slice of the simulation box parallel to the page. (c) sv20 droplet (in red) is depicted togetherwith the real part of solvent (dipole) polarization operator (cid:101) m ( r ) as arrows. (d)–(f): Same as (a)–(c), butat T ∗ = 50. Here, the polyampholytes and solvent are seen to distribute in an essentially uniform mannerthroughout the simulation box. The Distinctive Solvent Dielectric Environments in the Dilute and CondensedPhases Likely Contribute to a Minor Enhancement of Polyampholyte LLPSPropensity.
As reported above, our molecular dynamics and analytical RPA theoriesindicate, at least for the models considered, that allowing the solvent-contributed effective (cid:15) r inside a condensed polyampholyte droplet to decrease by using a physically plausible explicitpolar solvent model (Fig. 4) leads in most cases to an enhancement of polyampholyte LLPSpropensity relative to that predicted by commonly utilized implicit-solvent formulations inwhich the polyampholytes interact electrostatically via (cid:15) r of the bulk solvent (Figs. 6, 11,and 12). However, contrary to earlier suggestions that this effect can be dramatic, theenhancement effect observed here is only minor or at most moderate. We may quantify theenhancement effect by the ratio of the critical temperature of the explicit-solvent model tothat of the implicit-bulk-solvent, constant (cid:15) r model, or the percentage increase of the formerover the latter. For the molecular dynamics results in Fig. 6, this percentage increase incritical temperature, [ T ∗ , cr (ii) − T ∗ , cr (iii)] /T ∗ , cr (iii), obtained from the T ∗ , cr values in Table 1is approximately −
2% for sequence sv15, 5% for sv20, and 11% for sv30. These changesare minor, but nonetheless exhibit an increasing trend with increasing blockiness of thesequence charge pattern as characterized by − SCD, from essentially no change in T ∗ , cr forsv15 ( − SCD = 4 .
35) to a minor increase of 11% for sv30 ( − SCD = 27 . T ∗ , cr of the explicit-solvent model (ii) is almost identicalto that of the explicit nonpolar solvent model with (cid:15) r = 1 (i), suggesting in this casea relatively strong enhancement of polyampholyte LLPS by the heterogeneous dielectricenvironment entailed by the explicit-solvent model.The corresponding RPA results in Fig. 12a also show minor to at most moderateincreases. The T ∗ , cr (cid:15) r values for the explicit dipole, constant permittivity, and empiricalpermittivity models in Fig. 12a are, respectively, 4 .
46, 3 .
37, and 3 .
24 for sv15, 7 .
10, 6 . .
49 for sv20, and 23 .
7, 21 .
3, and 19 . T ∗ , cr (cid:15) r of the explicit dipole model over that of the constant permittivity model (differencebetween the two T ∗ , cr (cid:15) r values) tends to increase with − SCD (1 .
09 for sv15, 0 .
94 for sv20,and 2 .
40 for sv30), the perecentage increase in T ∗ , cr decreases with − SCD: approximately32% for sv15, 15% for sv20, and 11% for sv30. To what extent are these differences betweenthe molecular dynamics and RPA models attributable to their different bulk-solvent (cid:15) r values (4.0 for molecular dynamics, 20.0 for RPA) and their substantially differentcondensed-phase polymer volume fractions (much lower for RPA than for moleculardynamics) remains to be investigated. These uncertainties notwithstanding, the results inFigs. 6, 11, and 12 suggest that when all interactions among polyampholytes and polarsolvents are appropriately taken into account, the enhancement of polyampholyte LLPS bythe concentration-dependent dielectric heterogeneity of the polar solvent is likely minor.At the very least, one can be certain that the physical effect is much more modest than0when part of the solvent-polymer interactions are neglected as in the case of the unphysical“no self energy” model depicted by the dotted red curve in Fig. 12b. T *0 r sv30sv20sv15 Fig. 14:
Effective relative permittivity of condensed polyampholyte droplets in dipole solvent. Moleculardynamics simulation results are shown for polyampholyte sequences sv15 (circles), sv20 (squares) and sv30(diamonds) for the same systems in Fig. 4. Overall effective (cid:15) r inside each droplet (data points connectedby solid lines) is computed by applying Eq. 16 to local and global dipole moments contributed from boththe polyampholytes and the dipole solvent molecules, wherein the local dipole contribution is averagedover the volume of the droplet delimited by an interval of width ∆ z in the z -direction of the simulationbox. Included for comparison are effective (cid:15) r for each droplet contributed by the dipole solvent alone (datapoints connected by dashed lines). The widths ∆ z/a of the droplets used for this calculation, listed byincreasing temperature (same T ∗ s as those in Fig. 4), are: sv15: 40, 45, 60, 80; sv20: 40, 45, 60, 85; andsv30: 45, 60, 70, 80. The Dielectric Environment Experienced by Client Molecules in a Con-densed Polyampholyte Droplet Can be Complex.
We have focused so far on the (cid:15) r contributed by the polar solvent and its concentration dependence (Figs. 4, 8, and 9)in order to gauge the solvent’s role in the effective LLPS-driving electrostatic interactionsamong the polyampholytes. In these formulations, the charges on the polyampholytes aretreated explicitly, not as a part of the solvent-contributed dielectric background. Oncethe polyampholyte droplet—as a model biomolecular condensate—is formed, however, thedielectric environment entailed by the droplet as a whole—including contributions fromboth the polar solvent and the polyampholyte chains—is of potential biophysical andbiochemical interest because, among its many effects, it would affect how other molecules,1sometimes referred to as clients, partition into and interact within the droplet. Taking afirst step to address this interesting question, we compute, using our molecular dynamicsmodel, the overall effective relative permittivities of droplets that account for contributionsfrom the polar solvent as well as the polyampholytes themselves (Fig. 14). In the context ofbiomolecular modeling, this overall effective (cid:15) r would approximate the relative permittivityexperienced by client molecules (test charges) inside a biomolecular condensate. Fig. 14shows that this overall droplet (cid:15) r increases dramatically with − SCD of the sequences.For the three sequences considered, the overall droplet (cid:15) r s far exceed the model solvent’seffective (cid:15) r of 4 .
0. For sv30, the overall droplet (cid:15) r at low temperatures is seen to exceedeven the ≈
80 value for water. These high values of overall droplet (cid:15) r is likely caused bythe highly charged nature of the model polyampholytes we studied. Intuitively, one mayexpect less dramatic dielectric effects in biomolecular condensates scaffolded by naturalIDP polyampholytes that are less charged. This and related questions regarding theexpected spatial heterogeneity of the dielectric environment within condensed dropletsremain to be further explored. CONCLUSION
In summary, our explicit-chain molecular dynamics as well as analytical theory sug-gest that the heterogeneous solvent dielectric environment entailed by polyampholyteLLPS—with less solvent inside the condensed phase than outside—entails in most casesa minor but nonetheless appreciable enhancement of LLPS propensity relative to thatpredicted by polyampholytes interacting in a uniform dielectric environment carryingthe bulk-solvent relative permittivity. While these observations are based on simplifiedmodels, they suggest that the common approach of modeling polyampholyte LLPS with auniform (cid:15) ≈
80 likely leads to a minor underestimation of LLPS propensity. The extentof the expected small error, however, has to be ascertained by futher investigation, usingsimulations with realistic water models if possible. Further advances building on theanalytical formulation developed here should also provide a complementary approach toaddress this basic question. Although concentration-dependent solvent permittivity likelyincurs only a minor effect on the stability of polyampholyte condensates, we find that theoverall dielectric environment contributed by both the polar solvent and the polyampholytescaffold inside a condensate can be highly sensitive to the sequence charge pattern ofthe polyampholyte. Much about the complexity of the internal dielectric environment ofbiomolecular condensates and its biophysical and biochemical ramifications remain to beinvestigated.
Acknowledgements.
We thank Julie Forman-Kay for helpful discussions. Financial2support for this work was provided by Canadian Institutes of Health Research grantNJT-155930 and Natural Sciences and Engineering Research Council of Canada Discoverygrant RGPIN-2018-04351. We are grateful for the computational resources provided byCompute/Calcul Canada.The authors declare no conflict of interest.3
References Banani, S. F.; Lee, H. O.; Hyman, A. A.; Rosen, M. K. Biomolecular condensates: organizersof cellular biochemistry.
Nat. Rev. Mol. Cell Biol. , , 285–298. Hyman, A. A.; Weber, C. A.; J¨ulicher, F. Liquid-liquid phase separation in biology.
Annu.Rev. Cell Dev. Biol. , , 39–58. Brangwynne, C. P.; Tompa, P.; Pappu, R. V. Polymer physics of intracellular phase transitions.
Nat. Phys. , , 899–904. Lin, Y.-H.; Forman-Kay, J. D.; Chan, H. S. Sequence-specific polyampholyte phase separationin membraneless organelles.
Phys. Rev. Lett. , 2016, , 178101. Alberti, S.; Gladfelter, A.; Mittag, T. Considerations and challenges in studying liquid-liquidphase separation and biomolecular condensates.
Cell , , 419–434. Brangwynne, C. P.; Eckmann, C. R.; Courson, D. S.; Rybarska, A.; Hoege, C.; Gharakhani, J.;J¨ulicher, F.; Hyman, A. A. Germline P granules are liquid droplets that localize by controlleddissolution/condensation.
Science , , 1729–1732. Li, P.; Banjade, S.; Cheng, H. C.; Kim, S.; Chen, B.; Guo, L.; Llaguno, M.;Hollingsworth, J. V.; King, D. S.; Banani, S. F.; Russ, P. S.; Jiang, Q.-X.; Nixon, B. T.;Rosen, M. K. Phase transitions in the assembly of multivalent signalling proteins.
Nature , , 336–340. 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.; Frantz, 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: lowcomplexity sequence domains form dynamic fibers within hydrogels.
Cell , , 753–767. 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. Phase transition of a disordered nuageprotein generates environmentally responsive membraneless organelles.
Mol. Cell , ,936–947. 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 assemblyand drives pathological fibrillization.
Cell , , 123–133. 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 , , 208–219. Bergeron-Sandoval, L.P.; Safaee, N.; Michnick, S. W. Mechanisms and consequences of macro-molecular phase separation.
Cell , , 1067–1079. Shin, Y.; Brangwynne, C.P. Liquid phase condensation in cell physiology and disease.
Science , , eaaf4382. 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. , , 420–435. Ladouceur, A.-M.; Parmar, B. S.; Biedzinski, S.; Wall, J.; Tope, S. G.; Cohn, D.; Kim,A.; Soubry, N.; Reyes-Lamothe, R.; Weber, S. C. Clusters of bacterial RNA polymerase arebiomolecular condensates that assemble through liquid–liquidphase separation.
Proc. Natl.Acad. Sci. U.S.A. , , 18540–18549. Chong, P. A.; Forman-Kay, J. D. Liquid-liquid phase separation in cellular signaling systems.
Curr. Opin. Struct. Biol. , , 180–186. Li, X.-H.; Chavali, P. L.; Pancsa, R.; Chavali, S.; Babu, M. M. Function and regulation ofphase-separated biological condensates.
Biochemistry , , 2452–2461. Chen, X.; Wu, X.; Wu, H.; Zhang, M. Phase separation at the synapse.
Nat. Neurosci. , , 301–310. Gomes, E., and Shorter, J. The molecular language of membraneless organelles.
J. Biol. Chem. , , 7115–7127. Laflamme, G.; Mekhail, K. Biomolecular condensates as arbiters of biochemical reactionsinside the nucleus.
Commun. Biol. , , 773. Schaffert, N.; Hossbach, M.; Heintzmann, R.; Achsel, T.; L¨uhrmann, R. RNAi knockdown ofhPrp31 leads to an accumulation of U4/U6 di-snRNPs in Cajal bodies.
EMBO J. , ,3000–3009. Nott, T. J.; Craggs, T. D.; Baldwin, A. J. Membraneless organelles can melt nucleic acidduplexes and act as biomolecular filters.
Nat. Chem. , , 569–575. Lu, Y.; Wu, T.; Gutman, O.; Lu, H.; Zhou, Q.; Henis, Y. I.; Luo, K. Phase separation of TAZcompartmentalizes the transcription machinery to promote gene expression.
Nat. Cell Biol. , , 453–464. Csizmok, V.; Follis, A. V.; Kriwacki, R. W.; Forman-Kay, J. D. Dynamic protein interactionnetworks and new structural paradigms in signaling.
Chem. Rev. , , 6424–6462. Monahan, Z.; Ryan, V. H.; Janke, A. M.; Burke, K. A.; Rhoads, S. N.; Zerze, G. H.; O’Meally,R.; Dignon, G. L.; Conicella, A. E.; Zheng, W.; Best, R. B.; Cole, R. N.; Mittal, J.; Shewmaker,F.; Fawzi N. L. Phosphorylation of the FUS low-complexity domain disrupts phase separation,aggregation, and toxicity.
EMBO J. , , 2951–2967. Elbaum-Garfinkle, S.; Kim, Y.; Szczepaniak, K. Chen, C. C.-H.; Eckmann, C. R.; Myong, S.;Clifford P. Brangwynne, C. P. The disordered P granule protein LAF-1 drives phase separationinto droplets with tunable viscosity and dynamics.
Proc. Natl. Acad. Sci. U.S.A. , ,7189–7194. van der Lee, B.; Buljan, M.; Lang, B.; Weatheritt, R. J.; Daughdrill, G. W.; Dunker, A. K.;Fuxreiter, M.; Gough, J.; Gsponer, J.; Jones, D. T.; Kim, P. M.; Kriwacki, R. W.; Oldfield,C. J.; Pappu, R. V.; Tompa, P.; Uversky, V. N.; Wright, P. E.; Babu, M. M. Classification ofintrinsically disordered regions and proteins. Chem. Rev. , , 6589—6631. Lin, Y.-H.; Song, J.; Forman-Kay, J. D.; Chan, H. S. Random-phase-approximation theoryfor sequence-dependent, biologically functional liquid-liquid phase separation of intrinsicallydisordered proteins.
J. Mol. Liq. , 2017, , 176–193. Zhou, H.-X.; Pang, X. Electrostatic interactions in protein structure, folding, binding, andcondensation.
Chem. Rev. , , 1691–1741. Vernon, R. M.; Chong, P. A.; Tsang, B.; Kim, T. H.; Bah, A.; Farber, P.; Lin, H.; Forman-Kay, J. D. Pi-Pi contacts are an overlooked protein feature relevant to phase separation. eLife , , e31486. Wang, J.; Choi, J. M.; Holehouse, A. S.; Lee, H. O.; Zhang, X.; Jahnel, M.;Maharana, S.;Lemaitre, R.; Pozniakovsky, A.; Drechsel, D.; Poser, I.; Pappu, R. V.; Alberti, S.; Hyman, A.A. A molecular grammar governing the driving forces for phase separation of prion-like RNAbinding proteins.
Cell , , 688–699. Vernon, R. M.; Forman-Kay, J. D. First-generation predictors of biological protein phaseseparation.
Curr. Opin. Struct. Biol. , , 88–96. Martin, E. W.; Holehouse, A. S.; Peran, I.; Farag, M.; Incicco, J. J.; Bremer, A.; Grace, C. R.;Soranno, A.; Pappu, R. V.; Mittag, T. Valence and patterning of aromatic residues determinethe phase behavior of prion-like domains.
Science , , 694–699. Das, S.: Lin, Y.-H.; Vernon, R. M.; Forman-Kay, J. D.; Chan, H. S. Comparative roles ofcharge, π , and hydrophobic interactions in sequence-dependent phase separation of intrinsicallydisordered proteins. Proc. Natl. Acad. Sci. U.S.A. , , 28795–28805. Muiznieks, L. D.; Sharpe, S.; Pom`es, R.; Keeley, F. W. Role of liquid–liquid phase separationin assembly of elastin and other extracellular matrix proteins.
J. Mol. Biol. , , 4741–4753. Lin, Y.-H.; Forman-Kay, J. D.; Chan, H. S. Theories for sequence-dependent phase behaviorsof biomolecular condensates.
Biochemistry , , 2499–2508. Cinar, H.; Cinar, S.; Chan, H. S.; Winter, R. Pressure-induced dissolution and reentrantformation of condensed, liquid-liquid phase-separated elastomeric α -elastin. Chem. Eur. J. , 8286–8291. Dignon, G. L.; Zheng, W.; Kim, Y. C.; Mittal, J. Temperature-controlled liquid-liquid phaseseparation of disordered proteins.
ACS Cent. Sci. , , 821–830. Cinar, H.; Fetahaj, Z.; Cinar, S.; Vernon, R. M.; Chan, H. S.; Winter, R. Temperature, hydro-static pressure, and osmolyte effects on liquid-liquid phase separation in protein condensates:Physical chemistry and biological implications.
Chem. Eur. J. , , 13049–13069. Bianchi, G.; Longhi, S.; Grandori, R.; Brocca, S. Relevance of electrostatic charges in com-pactness, aggregation, and phase separation of intrinsically disordered proteins.
Int. J. Mol.Sci. , , 6208. Yeong, V.; Werth, E. G.; Brown, L. M.; Obermeyer, A. C. Formation of biomolecular conden-sates in bacteria by tuning protein electrostatics.
ACS Cent. Sci. , , 2301–2310. Tsang, B.; Arsenault, J.; Vernon, R. M.; Lin, H.; Sonenberg, N.; Wang, L. Y.; Bah, A.;Forman-Kay, J. D. Phosphoregulated FMRP phase separation models activity-dependenttranslation through bidirectional control of mRNA granule formation.
Proc. Natl. Acad. Sci.U.S.A. , , 4218–4227. Schuster, B. S.; Dignon, G. L.; Tang, W. S.; Kelley, F. M.; Ranganath, A. K.; Jahnke, C.N.; Simplins, A. G.; Regy, R. M.; Hammer, D. A.; Good, M. C.; Mittal, J. Identifying se-quence perturbations to an intrinsically disordered protein that determine its phase separationbehavior.
Proc. Natl. Acad. Sci. U.S.A. , , 11421–11431. Das, R. K.; Pappu, R. V. Conformations of intrinsically disordered proteins are influencedby linear sequence distributions of oppositely charged residues.
Proc. Natl. Acad. Sci. U.S.A. , , 13392–13397. Sawle, L.; Ghosh, K. A theoretical method to compute sequence dependent configurationalproperties in charged polymers and proteins.
J. Chem. Phys. , , 085101. Lin, Y.-H.; Chan, H. S. Phase separation and single-chain compactness of charged disorderedproteins are strongly correlated.
Biophys, J. , , 2043–2046. 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 propertiesof an intrinsically disordered region of a germ cell-specific protein on phase separation.
Proc.Natl. Acad. Sci. U.S.A. , , E8194–E8203. Amin, A. N.; Lin, Y.-H.; Das, S.; Chan, H. S. Analytical theory for sequence-specific binaryfuzzy complexes of charged intrinsically disordered proteins.
J. Phys. Chem. B , ,6709–6720. 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 disorderedproteins.
New J. Phys. , , 115003. Chang, L.-W.; Lytle, T. K.; Radhakrishna, M.; Madinya, J. J.; V´elez, J.; Sing, C. E.; Perry,S. L. Sequence and entropy-based control of complex coacervates.
Nat. Comm. , , 1273. Lin, Y.-H.; Brady, J. P.; Chan, H. S.; Ghosh, K. A unified analytical theory of heteropolymersfor sequence-specific phase behaviors of polyelectrolytes and polyampholytes.
J. Chem. Phys. , , 045102. Dignon, G. L.; Zheng, W.; Kim, Y. C.; Best, R. B.; Mittal, J. Sequence determinants ofprotein phase behavior from a coarse-grained model.
PLoS Comput. Biol. , , e1005941. Das, S.; Eisen, A.; Lin, Y.-H.; Chan, H. S. A lattice model of charge-pattern-dependentpolyampholyte phase separation.
J. Phys. Chem. B , , 5418–5431. 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. , , 9929–9934. Das, S.; Amin, A. N.; Lin, Y.-H.; Chan, H. S. Coarse-grained residue-based models of disor-dered protein condensates: Utility and limitations of simple charge pattern parameters.
Phys.Chem. Chem. Phys. , , 28558–28574. Hazra, M. K.; Levy, Y. Charge pattern affects the structure and dynamics of polyampholytecondensates.
Phys. Chem. Chem. Phys. , , 19368–19375. Robichaud,N. A. S.; Saika-Voivod, I.; Wallin, S. Phase behavior of blocky charge lattice poly-mers: Crystals, liquids, sheets, filaments, and clusters.
Phys. Rev. E. , , 052404. 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 undercellular conditions eLife , , e42571. McCarty, J.; Delaney, K. T.; Danielsen, S. P. O.; Fredrickson, G. H.; Shea, J.-E. Completephase diagram for liquid-liquid phase separation of intrinsically disordered proteins.
J. Phys.Chem. Lett. , , 1644–1652. Danielsen, S. P. O.; McCarty, J.; Shea, J.-E.; Delaney, K. T.; Fredrickson, G. H. Moleculardesign of self-coacervation phenomena in block polyampholytes.
Proc. Natl. Acad. Sci. U.S.A. , , 8224–8232. Danielsen, S. P. O.; McCarty, J.; Shea, J.-E.; Delaney, K. T.; Fredrickson, G. H. Small ioneffects on self-coacervation phenomena in block polyampholytes.
J. Chem. Phys. , ,034904. Pal, T.; Wess´en, J.; Das, S.; Chan, H. S. Subcompartmentalization of polyampholyte speciesin organelle-like condensates is promoted by charge pattern mismatch and strong excluded-volume interaction. , available on arXiv: https://arxiv.org/abs/2006.12776 . Gilson, M. K.; Honig, B. H. The dielectric constant of a folded protein.
Biopolymers , , 2097–2119. Garc´ıa-Moreno, B.; Dwyer, J.; Gittis, A. G.; Lattman, E. E.; Spencer, D. S.; Stites, W. E.Experimental measurement of the effective dielectric in the hydrophobic core of a protein.
Biophys. Chem. , , 211–224. Dwyer, J. J.; Gittis, A. G.; Karp, D. A.; Lattman, E. E.; Spencer, D. S.; Stites, W. E.;Garc´ıa-Moreno E, B. High apparent dielectric constants in the interior of a protein reflectwater penetration.
Biophys. J. , , 1610–1620. Pitera, J. W.; Falta, M.; van Gunsteren, W. F. Dielectric properties of proteins from simula-tion: The effects of solvent, ligands, pH, and temperature.
Biophys. J . , , 2546-2555. Schutz, C. N.; Warshel, A. What are the dielectric “constants” of proteins and how to validateelectrostatic models?
Proteins , , 400–417. Guest, W. C.; Cashman, N. R.; Plotkin, S. S. A theory for the anisotropic and inhomogeneousdielectric properties of proteins.
Phys. Chem. Chem. Phys. , , 6286–6295. Cuervo, A.; Dans, P. D.; Carrascosa, J. L.; Orozco, M.; Gomila, G.; Fumagalli, L. Directmeasurement of the dielectric polarization properties of DNA.
Proc. Natl. Acad. Sci. U.S.A. , , E3624–E3630. Fumagalli, L.; Esfandiar, A.; Fabregas, R.; Hu1, S.; Ares, P.; Janardanan, A.; Yang, Q.;Radha1, B.; Taniguchi, T.; Watanabe, K. Gomila, G.; Novoselov, K. S.; Geim, A. K. Anoma-lously low dielectric constant of confined water.
Science , , 1339–1342. Bragg, W. L.; Pippard, A. B. The form birefringence of macromolecules.
Acta Cryst . , ,865—867. Jackson, J. D.
Classical Electrodynamics (2nd edition); Wiley, :New York, 1975; pp 154–155. Markel, V. A. Introduction to the Maxwell Garnett approximation: Tutorial.
J. Opt. Soc. Am.A Ferreira, L. A.; Loureiro, J. A.; Gomes, J.; Uversky, V. N.; Madeira, P. P.; Zaslavsky, B.Y. Why physicochemical properties of aqueous solutions of various compounds are linearlyinterrelated.
J. Mol. Liq. , , 116–123. Moreau, N. G.; Martin, N.; Gobbo, P.; Tang, T.-Y. D.; Mann, S. Spontaneous membrane-lessmulti-compartmentalization via aqueous two-phase separation in complex coacervate micro-droplets.
Chem. Commun. , , 12717. Yildirim, A.; Sharma, M.; Varner, B. M.; Fang, L.; Feig, M. Conformational preferences ofDNA in reduced dielectric environments.
J. Phys. Chem. B , , 10874–10881. Arscott, P. G.; Ma, C.; Wenner, J. R.; Bloomfield, V. A. DNA condensation by cobalt hexaam-mine(III) in alcohol-water mixtures: Dielectric constant and other solvent effects.
Biopolymers , , , 345–364. Nakashima, K. K.; Vibhute, M. A.; Spruijt, E. Biomolecular chemistry in liquid phase sepa-rated compartments.
Frontiers Mol. Biosci. , , 21. Milin, A. N.; Deniz, A. A. Reentrant phase transitions and non-equilibrium dynamics inmembraneless organelles.
Biochemistry , , 2470–2477. Zaslavsky, B. Y.; Uversky, V. N. In aqua veritas: The indispensable yet mostly ignored role ofwater in phase separation and membrane-less organelles.
Biochemistry , , 2437–2451. Zaslavsky, B. Y.; Ferreira, L. A.; Darling, A. L.; Uversky, V. N. The solvent side of proteina-ceous membrane-less organelles in light of aqueous two-phase systems.
Int. J. Biol. Macromol. , , 1224–1251. Pliss, A.; Levchenko, S. M.; Liu, L.; Peng, X.; Ohulchanskyy, T. Y.; Roy, I.; Kuzmin, A.N.; Qu, J.; Prasad, P. N. Cycles of protein condensation and discharge in nuclear organelles studied by fluorescence lifetime imaging. Nat. Comm. , , 455. Zheng, W.; Dignon, G. L.; Jovic, N.; Xu, X.; Regy, R. M.; Fawzi, N. L.; Kim, Y. C.; Best, R.B.; Mittal, J. Molecular details of protein condensates probed by microsecond long atomisticsimulations.
J. Phys. Chem. B , , 11671–11679. Benayad, Z.; von B¨ulow, S.; Stelzl, L. S.; Hummer, G. Simulation of FUS protein con-densates with an adapted coarse-grained model.
J. Chem. Theory Comput. , doi:10.1021/acs.jctc.0c01064 de Jong, D. H.; Singh, G.; Bennett, W. F. D.; Arnarez, C.; Wassenaar, T. A.; Sch¨afer, L.V.; Periole, X.; Tieleman, D. P.; Marrink, S. J. Improved parameters for the Martini coarse-grained protein force field. J. Chem. Theory Comput. , , 687–697. Michalowsky, J.; Sch¨afer, L. V.; Holm, C.; Smiatek, J. A refined polarizable water model forthe coarse-grained MARTINI force field with long-range electrostatic interactions.
J. Chem.Phys . , , 054501. K¨ofinger, J.; Hummer, G.; Dellago, C. A one-dimensional dipole lattice model for water innarrow nanopores.
J. Chem. Phys. , , 154110. Motevaselian, M. H.; Mashayak, S. Y.; R. Aluru, N. R. Extended coarse-grained dipole modelfor polar liquids: Application to bulk and confined water.
Phys. Rev. E , , 052135. Levy, A.; Andelman, D.; Orland, H. Dielectric constant of ionic solutions: A field-theoryapproach.
Phys. Rev. Lett. , , 227801. Levy, A.; Andelman, D.; Orland, H. Dipolar Poisson-Boltzmann approach to ionic solutions:A mean field and loop expansion analysis.
J. Chem. Phys. , , 164909. Zhuang, B.; Wang, Z.-G. Statistical field theory for polar fluids.
J. Chem. Phys. , , Martin, J. M.; Li, W.; Delaney, K. T.; Fredrickson, G. H. Statistical field theory descriptionof inhomogeneous polarizable soft matter.
J. Chem. Phys. , , 154104. Grzetic, D. J.; Delaney, K. T.; Fredrickson, G. H. Constrasting dielectric properties of elec-trolyte solutions with polar and polarizable solvents.
Phys. Rev. Lett. , , 128007. Martin, J. M.; Delaney, J. T.; Fredrickson, G. H. Effect of an electric field on the stability ofbinary dielectric fluid mixtures.
J. Chem. Phys. , , 234901. Silmore, K. S.; Howard, M. P.; Panagiotopoulos, A. Z. Vapor-liquid equilibrium and surfacetension of fully flexible Lennard-Jones chains.
Mol. Phys. , , 320–327. Mundy, C. J.; Siepmann, J. I.; Klein, M. L. Calculation of the shear viscosity of decane usinga reversible multiple timestep algorithm.
J. Chem. Phys. , , 3376–3380. Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. United-atomdescription of n-alkanes.
J. Phys. Chem. B , , 2569–2577. Nicolas, J. P.; Smit, B. Molecular dynamics simulations of the surface tension of n-hexane,n-decane and n-hexadecane.
Mol. Phys. , , 2471–2475. P`amies, J. C.; McCabe, C.; Cummings, P. T.; Vega, L. F. Coexistence densities of methane and propane by canonical molecular dynamics and gibbs ensemble Monte Carlo simulations. Mol. Simul. , , 463–470. Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining theequilibrium structure of simple liquids.
J. Chem. Phys. , , 5237–5247. Anderson, J. A.; Lorenz, C. D.; Travesset, A. General purpose molecular dynamics simulationsfully implemented on graphics processing units.
J. Comput. Phys. , , 5342–5359. Glaser, J.; Nguyen, T. D.; Anderson, J. A.; Lui, P.; Spiga, F.; Millan, J. A.; Morse, D. C.;Glotzer, S. C. Strong scaling of general-purpose molecular dynamics simulations on GPUs.
Comput. Phys. Comm. , , 97–107. LeBard, D. N.; Levine, B. G.; Mertmann, P.; Barr, S. A.; Jusufi, A.; Sanders, S.; Klein, M.L.; Panagiotopoulos, A. Z. Self-assembly of coarse-grained ionic surfactants accelerated bygraphics processing units.
Soft Matter , , 2385–2397. Mart´ınez, L.; Andrade, R.; Birgin, E. G.; Mart´ınez, J. M. PACKMOL: A package for buildinginitial configurations for molecular dynamics simulations.
J. Comput. Chem. , , 2157–2164. Abrashkin, A.; Andelman, D.; Orland, H. Dipolar Poisson-Boltzmann equation: Ions anddipoles close to charged interfaces.
Phys. Rev. Lett. , , 077801. Ballenegger, V.; Hansen, J.-P. Dielectric permittivity profiles of confined polar fluids.
J. Chem.Phys. , , 114711. Stern, H. A.; Feller, S. E. Calculation of the dielectric permittivity profile for a nonuniformsystem: Application to a lipid bilayer simulation.
Zhu, H.; Yang, F.; Zhu, Y.; Li, A.; He, W.; Huang, J.; Li, G. Investigation of dielectricconstants of water in a nano-confined pore.
RSC Adv. , , , 8628–8635. J. Chem. Phys. , , 3401–3412. Chan, H. S.; Dill, K. A. Solvation: How to obtain microscopic energies from partitioning andsolvation experiments.
Annu. Rev. Biophys. Biomol. Struct. , , 425–459. Ray, A. Solvophobic interactions and micelle formation in structure forming nonaqueous sol-vents.
Nature , , 313–315. Wang, Z.-G. Fluctuation in electrolyte solutions: The self energy.
Phys. Rev. E , ,021501. Chandler, D. The dielectric constant and related equilibrium properties of molecular fluids:Interaction site cluster theory analysis.
J. Chem. Phys. , , 1113–1124. Edwards, S. F. The statistical mechanics of polymers with excluded volume.
Proc. Phys. Soc. , , 613–624. Fredrickson, G. H.
The Equilibrium Theory Of Inhomogeneous Polymers ; Oxford UniversityPress Inc.:New York, 2006.
Fredrickson, G. H.; Ganesan, V.; Drolet, F. Field-theoretic computer simulation methods for polymers and complex fluids. Macromolecules , , 16–39. Parisi, G. On complex probabilities.
Phys. Lett. B , , 393–395. Klauder, J. R. A Langevin approach to fermion and quantum spin correlation functions.
J.Phys. A: Math. Gen. , , L317–L319. Parisi, G.; Wu, Y.-S. Perturbation theory without gauge fixing.
Scientia Sinica ,483–496. Chan, H. S.; Halpern, M. B. New ghost-free infrared-soft gauges.
Phys. Rev. D , ,540–547. Chan, H. S.; Halpern, M. B. Continuum-regularized quantum gravity.
Zeitschrift F¨ur PhysikC , , 669–693. Lennon, E. M.; Mohler, G. O.; Ceniceros, H. D.; Garca-Cervera, C. J.; Fredrickson, G. H.Numerical solutions of the complex Langevin equations in polymer field theory.
MultiscaleModeling & Simulation , , 1347–1370. Delaney, K. T.; Fredrickson, G. H. Theory of polyelectrolyte complexation–complex coacer-vates are self-coacervates.
J. Chem. Phys. , , 224902. Riggleman, R. A.; Fredrickson, G. H. Field-theoretic simulations in the gibbs ensemble.
J.Chem. Phys. , , 024104.2