Decoding the physical principles of twocomponent biomolecular phase separation
Abstract
Cells possess a multiplicity of nonmembranebound compartments, which form via liquidliquid phase separation. These condensates assemble and dissolve as needed to enable central cellular functions. One important class of condensates is those composed of two associating polymer species that form onetoone specific bonds. What are the physical principles that underlie phase separation in such systems? To address this question, we employed coarsegrained molecular dynamics simulations to examine how the phase boundaries depend on polymer valence, stoichiometry, and binding strength. We discovered a striking phenomenon – for sufficiently strong binding, phase separation is suppressed at rational polymer stoichiometries, which we termed the magicratio effect. We further developed an analytical dimergel theory that confirmed the magicratio effect and disentangled the individual roles of polymer properties in shaping the phase diagram. Our work provides new insights into the factors controlling the phase diagrams of biomolecular condensates, with implications for natural and synthetic systems.
Introduction
Eukaryotic cells are host to a multiplicity of nonmembranebound compartments. Recent studies have shown that these compartments form via liquidliquid phase separation (Brangwynne et al., 2009; Li et al., 2012; Molliex et al., 2015). The phaseseparated condensates enable many central cellular functions – from ribosome assembly, to RNA regulation and storage, to signaling and metabolism (Shin and Brangwynne, 2017; Banani et al., 2017). Unlike conventional liquidliquid phase separation, for example wateroil demixing, the underlying interactions that drive biomolecular phase separation typically involve strong onetoone saturable interactions, often among multiple components (Ditlev et al., 2018). As a result, the phase diagrams of biomolecular condensates are complex and are sensitive to a variety of physical properties of the biomolecules, included number of binding sites, binding strengths, and additional nonspecific interactions. Importantly, these physical parameters can be subject to biological regulation, and can thus directly impact the organization and function of the condensates. It is therefore crucial to understand how the physical properties of the components shape the phase diagram of biomolecular condensates.
Biomolecular condensates typically contain tens to hundreds of types of molecules. Yet, when characterized in detail, only a small number of components are responsible for condensate formation (Ditlev et al., 2018). One class of such condensates are those formed by the association of two essential components. In the simplest case, each component consists of repeated domains/stickers that bind in a onetoone fashion with the domains of the other component (Figure 1A and B; Choi et al., 2019; Xu et al., 2020). Such twocomponent condensates have been observed in both natural and engineered contexts. For example, the pyrenoid, an organelle responsible for carbon fixation in the alga Chlamydomonas reinhardtii, is a condensate of the CO_{2}fixing enzyme Rubisco with the linker protein Essential PYrenoid Component 1 (EPYC1). EPYC1 consists of five evenlyspaced Rubiscobinding regions, while Rubisco holoenzyme has eight specific binding sites for EPYC1. Multivalent interactions between Rubisco and EPYC1 are responsible for pyrenoid formation (Freeman Rosenzweig et al., 2017; Wunder et al., 2018; He et al., 2020). Promyelocytic leukemia (PML) nuclear bodies are condensates of PML proteins. PML is SUMOylated at three main positions and several minor sites. These modifications and a Cterminal SUMO Interaction Motif (SIM) found in most PML isoforms contribute to the formation of these bodies (Shen et al., 2006). Engineered polySUMO and polySIM proteins (10 repeats of Small Ubiquitinlike Modifier [SUMO] and SIM, respectively) phase separate when mixed together, but not as individual components (Banani et al., 2016; Ditlev et al., 2018).
Previous simulations (Freeman Rosenzweig et al., 2017; Xu et al., 2020) of average cluster size in such twocomponent systems revealed a striking phenomenon – for sufficiently strong binding, the formation of large clusters is suppressed when the valence of one species equals or is an integral multiple of the valence of the other species, favoring the formation of small stable oligomers instead of a condensate. The phenomenon reminiscent of the exact filling of atomic shells leading to the unreactive noble gases was termed the ‘magicnumber’ effect. A similar effect was found in a ternary system modeling the clustering of nephrin, Nck, and NWASP proteins which regulates cellcell adhesion in podocyte cells of the kidney (Chattaraj et al., 2019). However, cluster size may reflect a solgel percolation transition rather than a thermodynamic phase transition (Harmon et al., 2017), and thus provides at best a qualitative measure of phase separation. Moreover, these previous studies focused on equal sticker stoichiometry, whereas biomolecular condensates cover a broad range of stoichiometries both in vitro (Li et al., 2012; Banani et al., 2016) and in vivo (Sanders et al., 2020).
Here, we directly delineate the full phase diagram of such twocomponent systems. Using coarsegrained molecular dynamics simulations, we explore systematically how phase boundaries depend on valence, stoichiometry, and binding strength of two associating polymers (Figure 1C and D). Our studies reveal an unanticipated effect – when the numbers of polymers of the two types have a rational stoichiometry (1:1, 1:2, etc.), phase separation can be strongly suppressed, which we call the ‘magicratio’ effect (Figure 1D, phase diagram at low temperatures). To understand the magicratio effects better, we develop a twocomponent sticker theory à la Semenov and Rubinstein (Semenov and Rubinstein, 1998). We model the system as dominated by polymer dimers in the dilute phase and by a condensate of independent stickers in the dense phase (Figure 1B). The resulting analytical theory captures the magicratio effect discovered in simulations, and allows us to disentangle the individual roles of valence, stoichiometry, specificbond strength, and nonspecific attraction in determining the phase boundaries of twocomponent multivalent systems. Living cells regulate the valence and interactions of biomolecules through chemical modification, or on a slower timescale, tune the stoichiometry via synthesis/degradation or sequestration, and over evolutionary time, adapt the strength of specific and nonspecific interactions through mutation of molecular sequences. Understanding the individual roles of these biologically tunable variables thus brings new insights into possible cellular strategies for regulating the formation and dissolution of biomolecular condensates.
Results
Coarsegrained moleculardynamics simulations
We perform coarsegrained moleculardynamics simulations using LAMMPS (Plimpton, 1995) to determine the phase boundaries of twocomponent multivalent systems (Figure 2). Briefly, we model the two polymer species as flexible linear chains of beads connected by harmonic springs (Figure 2A). Each bead represents one associative domain/sticker of the polymer. To ensure associative domains of different polymer types bind in a onetoone fashion, we impose a finiteranged attractive interaction between beads of different types. This, however, could lead to more than onetoone associations. Therefore, to avoid such unwanted associations, we impose strong repulsive interactions between beads of the same type over a large enough range to prevent other beads overlapping with a bound pair, thus preventing multipletoone binding (Figure 2B and Appendix 1—figure 1), see Appendix 1 for details.
To find the binodal phase boundaries, we simulate hundreds of polymers of types A and B with, respectively, m and n stickers (an $\text{A}}_{m}:{\text{B}}_{n$ system) in a box with periodic boundary conditions (Figure 2C and D). We initialize the system by constructing a dense slab of polymers in the middle of the box (Dignon et al., 2018). The system evolves and relaxes according to Langevin dynamics (Langevin, 1908). After the system has achieved equilibrium, two phases coexist: a dilute phase consisting of dimers and other small oligomers, and a dense phase of an interconnected polymer condensate. We measure the corresponding density profile (Figure 2E) and calculate the dilute and densephase concentrations by averaging the density profile over the regions ($x\le 100\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$ or $x\ge 100\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$) and ($10\text{nm}\le x\le 10\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$), respectively. See Appendix 1 for simulation details.
Effect of valence
It was shown previously that for equal sticker stoichiometry in the strongbinding regime, clustering is substantially suppressed when the number of binding sites on one polymer species is an integer multiple of the number of binding sites on the other, as this condition favors the assembly of small oligomers in which all binding sites are saturated (Freeman Rosenzweig et al., 2017; Xu et al., 2020). What does this magicnumber effect imply for the actual phase diagram? To address this question, we fix the valence of polymer A at 14 and systematically vary the valence of polymer B from 5 to 16 while keeping the two sticker concentrations the same, that is, at equal global sticker stoichiometry.
Figure 3A and B show simulation results for the total sticker concentrations of the dilute and dense phases for A_{14}:B_{5} to A_{14}:B_{16} systems. In the strong binding regime, for magicnumber cases, that is when the valence of B is 7 or 14, the dilutephase concentration shows pronounced peaks (Figure 3A, black curve). What is the origin of the peak at A_{14}:B_{14}? Intuitively, when the dilute phase of the twocomponent system is dominated by dimers (for systems A_{14}:B_{12} to A_{14}:B_{16}, as supported by cluster size analysis in Appendix 1—figure 2), each of these dimers has high translational entropy, whereas polymers in the dense condensate have low translational entropy. For A_{14}:B_{14}, all binding sites can pair up in a dimer just as well as in the condensate, so the energy per polymer is not necessarily lower in the condensate. Why then is the condensate still competitive with the dilute phase? In a dimer, the binding sites of A_{14} must match all the binding sites of B_{14}, leading to a reduced overall conformational entropy. By comparison, the polymers in the condensate are more independent, binding to multiple members of the other species and enjoying a relatively higher overall conformational entropy. Because the translational entropy of each dimer decreases as their concentration goes up, the condensed phase eventually becomes more favorable and so the system phase separates with increasing concentration. Therefore, phase separation in A_{14}:B_{14} is primarily driven by a competition between translational entropy and conformational entropy.
By contrast, for A_{14}:B_{13} and A_{14}:B_{15}, one of the stickers in the dimer cannot be paired, and for A_{14}:B_{12} and A_{14}:B_{16}, two stickers per dimer cannot be paired. Therefore, forming a condensate not only increases the conformational entropy but more importantly lowers the energy of these systems. This significantly tilts the balance in favor of condensation. As a result, the dilutephase concentration is sharply peaked at A_{14}:B_{14}, falling off rapidly for increasingly unequal polymer lengths. We note that the densephase concentration shows no such feature (Figure 3B), indicating that the peak at A_{14}:B_{14} does not arise from differences in the internal structure of the dense phase.
The dilute phase of twocomponent systems is not always dominated by dimers (Appendix 1—figure 2). For example, the dilute phase of the A_{14}:B_{7} system is dominated by fullybonded trimers with 1 A_{14} and 2 B_{7}, the dilute phase of A_{14}:B_{8} is dominated by trimers with 1 A_{14} and 2 B_{8}, which has two unpaired stickers per trimer, and the dilute phase of A_{14}:B_{6} is dominated by oligomers with 3 A_{14} and 7 B_{6}, which although fullybonded is not small (Figure 3D). Consistent with the above logic, we find another peak in the dilutephase concentration at A_{14}:B_{7} (Figure 3A). More generally, in contrast to the magicnumber systems, the dilute phases in other cases are dominated by oligomers which are not capable of being fully bonded (high energy) and/or not small (low translational entropy) (Appendix 1—figure 2). The dilutephase concentration is therefore lower in these nonmagicnumber cases.
Effect of binding strength
How do the phase boundaries depend on the strength of binding? Figure 3A shows that, for nonmagicnumber systems, the dilutephase concentration decreases monotonically with increasing binding strength, whereas for magicnumber systems the dependence can be nonmonotonic. This difference is attributed to the distinct underlying driving forces for phase separation. For nonmagicnumber systems, as clustering allows a larger fraction of binding sites to be paired, the stronger the binding, the more the energy is lowered by condensate formation. Therefore, the dilutephase concentration drops as binding strength increases (or as temperature decreases). Such energydependence is expected for conventional phaseseparation models, such as Flory, 1942; Huggins, 1941.
Interestingly, for the magicnumber system A_{14}:B_{14}, the dilutephase concentration first decreases with increasing binding strength in the weak binding regime, similar to nonmagicnumber systems. However, as the binding energy is increased further, most of binding sites pair up in both dilute and dense phases. Phase separation is then primarily driven by a competition between conformational and translational entropy. The pairing up of binding sites reduces the conformational entropy of both the dense and dilute phases. By contrast, the translational entropy of the dilutephase components is almost unaffected. Consequently, the dilute phase becomes more competitive relative to the condensate, so the dilute phase boundary shifts to higher concentration.
By comparison, the densephase concentration increases monotonically with increasing binding strength for all systems (Figure 3B). This follows because the stronger the binding, the more stickers are paired, which tightens the condensate structure. We note that at substantially higher binding energies than studied here, essentially all the binding sites are satisfied in both magicnumber and nonmagicnumber systems, and the phase boundaries become independent of binding energy.
Effect of sticker stoichiometry
How do the phase boundaries depend on overall sticker stoichiometry? Figure 4A and B show total sticker concentrations of the dilute and dense phases for magicnumber systems A_{8}:B_{8} to A_{14}:B_{14} at different global sticker stoichiometries. For each system, the dilutephase concentration peaks at equal sticker ratio, falls off initially as the ratio deviates from 1, and then curves back up. What is the origin of the peak at equal sticker stoichiometry? Recall that, in the strong binding regime, phase separation of magicnumber systems is primarily driven by a competition between translational entropy and conformational entropy. Now consider starting with a system at equal sticker concentration, and adding more of one polymer species to the system. At the beginning, the added polymers readily enter the dense phase, which relaxes the conformational constraint that every sticker in the condensate has to pair with a partner. This increase of the conformational entropy of the condensate makes it more competitive, so the dilutephase concentration decreases. However, as the ratio between the two polymers is increased further, it becomes possible to form a spectrum of dilutephase oligomers which typically contain one extra polymer of the majority type (Appendix 1—table 1). These new oligomers have more relaxed structures than fully bonded dimers, which raises the conformational entropy of the dilute phase. Therefore, the dilute phase is favored over the condensate and its concentration curves back up.
Figure 4A also reveals that the dilutephase concentration decreases with increasing polymer valence. This follows in part because translational entropy in the dilute phase is per dimer center of mass, whereas conformational entropy in both phases scales with the number of stickers. The entropic gain of joining the dense phase is therefore more on a per sticker basis for longer polymers, so the dilutephase concentration decreases with increasing valence. As a less apparent yet important point, Figure 4A also shows that increasing polymer valence enhances both the width and relative height of the peak in the dilutephase concentration. The inferred phase diagram for the A_{8}:B_{8} system at ${U}_{0}=14{k}_{\text{B}}T$ is shown in Appendix 1—figure 3 together with the homogeneous gelation/percolation threshold obtained at ${U}_{0}=8{k}_{\text{B}}T$. We also report in Appendix 1—figure 5C the volume fraction of the polymers in the dense phase, which is ∼10%, comparable to the volume fraction of proteins in the cell cytoplasm.
Figure 4C and D show total sticker concentrations of the dilute and dense phases for unequal valence polymers A_{8}:B_{7} to A_{14}:B_{13} at different global sticker stoichiometries. The dilute phase boundary shows a symmetric minimum around equal stoichiometry for A_{8}:B_{7}, yet surprisingly, the phase boundary becomes asymmetric and then peaks at equal polymer stoichiometry with increasing polymer length (Figure 4C). What is the origin of these peaks? Taking the A_{14}:B_{13} system as an example, its dilute phase is dominated by dimers with an unpaired A sticker. This strongly disfavors the dilute phase in the strong binding regime at equal sticker stoichiometry. However, as the overall A:B sticker stoichiometry increases, the excess As cannot be paired anyway. In particular, at equal polymer stoichiometry (denoted as black dots in Figure 4C), forming dimers is no longer energetically costly. Therefore, to the left of the A_{14}:B_{13} peak at equal polymer stoichiometry, the dilutephase concentration is low because dimers are energetically disfavored as more bonds can be satisfied in the condensate. By contrast, to the right of the peak, the dilutephase concentration is low for a different reason – because the condensate is entropically favored, similar to the peak with respect to stoichiometry for magicnumber systems. Eventually, the dilutephase concentration curves back up due to formation of higher oligomers in the dilute phase, as discussed for magicnumber systems.
We note that for all these systems the densephase concentration shows no such striking features. Rather, the concentration decreases monotonically as the global sticker stoichiometry departs from one and as the valence of polymers decreases (Figure 4B and D).
Effect of valence and stoichiometry
Above, we considered the role of both relative valence and relative stoichiometry. By plotting phase boundaries as joint functions of valence and stoichiometry, we obtain a unified picture: Figure 5A and B show the dilute and densephase concentrations for systems A_{14}:B_{1216} at global sticker stoichiometries 14:1216. Notably, the dilutephase concentration is peaked along the diagonal (Figure 5A), that is at equal polymer stoichiometry, which we term the 'magicratio' effect because it occurs for rational ratios of associative polymers. Intuitively, all cases along the diagonal favor 1:1 polymer dimers: the dimers enjoy high translation entropy and there is no energy penalty involved in their formation. Thus, a dilute phase of dimers is strongly favored at equal polymer stoichiometry.
As for the dense phase concentration, it decreases monotonically as the global sticker stoichiometry departs from one and as the valence of polymers decreases (Figure 5B). This again indicates that the anomalous dependence of the dilutephase concentration on valence and stoichiometry does not arise from special properties of the dense phase.
Dimergel theory
While our simulations have revealed that a magicratio effect influences the boundaries of phase separation for associating polymers, we desire a deeper understanding of the interplay of factors such as overall valence, stoichiometry, and interaction strength. To this end, we develop a mean field theory of twocomponent associative polymers à la Semenov and Rubinstein (Semenov and Rubinstein, 1998; Xu, 2018).
Specifically, we consider a system of A and B polymers as in our simulations. Each polymer is a linear chain of L_{1} or L_{2} stickers of type A or type B, respectively. Without loss of generality, we take ${L}_{1}\ge {L}_{2}$. stickers of different types associate in a onetoone fashion. Our simulations suggest that for polymers of similar valence close to equal polymer stoichiometry the dilute phase is dominated by dimers and the dense phase is a gel network. Therefore, we assume that polymers can associate either as dimers or, alternatively, as a condensate in which pairs of stickers bind independently. This assumption of independence is a mean field approximation, as it neglects correlations between stickers in the same chain, and thus only applies when the polymers strongly overlap, that is at densities above the semidilute regime (De Gennes, 1979).
The partition function of such a system can be divided into three parts: $Z={Z}_{\text{ni}}{Z}_{\text{s}}{Z}_{\text{ns}}$, where ${Z}_{\text{ni}}$, the partition function of a solution of noninteracting polymers, captures the translational and conformational entropy of the two polymer species, ${Z}_{\text{s}}$ captures specific interactions between associating stickers, and ${Z}_{\text{ns}}$ captures all nonspecific interactions.
The corresponding freeenergy density for the mixed noninteracting polymers is Semenov and Rubinstein, 1998:
where c_{1} and c_{2} are the concentrations of A and B polymers measured in terms of stickers. Note that the terms for the conformational entropy of noninteracting polymers are omitted in Equation 1, as they are linear in c_{1} and c_{2} and thus do not influence the phase boundaries.
To include specific interactions, we first consider the partition function ${Z}_{\text{s}}({N}_{\text{d1}},{N}_{\text{d2}},{N}_{\text{b}})$ for states with exactly $N}_{\text{d1}$ and ${N}_{\text{d2}}$ total numbers of stickers of A and B types in dimers (i.e. number of dimers equals ${N}_{\text{d1}}/{L}_{1}={N}_{\text{d2}}/{L}_{2}$) and ${N}_{\text{b}}$ additional sticker pairs,
In Equation 2, P is the number of different ways that polymers and stickers can be chosen to pair up to form dimers and independent bonds,
where N_{1} and N_{2} are the total numbers of stickers of A and B types. (Note that in Equation (3) if ${L}_{1}>{L}_{2}$, the excess stickers of type A in dimers do not form additional bonds.) In Equation (2), W is the probability that all chosen polymers and stickers are, respectively, close enough to their specified partners in the noninteracting state to form dimers and independent bonds,
where ${v}_{\text{d}}$ and ${v}_{\text{b}}$ are effective interaction volumes and V is the system volume. The last term in Equation (2) is the Boltzmann factor for specific interactions, where ${\u03f5}_{\text{d}}$ and ${\u03f5}_{\text{b}}$ are the effective binding energies of dimers and sticker pairs, in units of ${k}_{\text{B}}T$.
The part of the freeenergy density due to specific interactions is
Using Stirling’s approximation $\mathrm{ln}N!=N\mathrm{ln}NN$, we obtain
where ${K}_{\text{d}}\equiv {e}^{{\u03f5}_{\text{d}}}/{v}_{\text{d}}$ and ${K}_{\text{b}}\equiv {e}^{{\u03f5}_{\text{b}}}/{v}_{\text{b}}$ are, respectively, the dissociation constants of a dimer and of a pair of stickers. ${c}_{\text{d1}}$ and ${c}_{\text{d2}}$ are the concentrations of stickers of A and B types in dimers (so ${c}_{\text{d1}}/{L}_{1}={c}_{\text{d2}}/{L}_{2}$), and ${c}_{\text{b}}$ is the concentration of independent bonds.
In the thermodynamic limit, ${F}_{\text{s}}$ will be minimized with respect to ${c}_{\text{d1}}$, ${c}_{\text{d2}}$ and ${c}_{\text{b}}$, which implies
Note that if ${c}_{\text{b}}$ in Equation (7) and ${c}_{\text{d1}}$ and ${c}_{\text{d2}}$ in Equation (8) are set to zero, these equations reduce to
where ${\rho}_{1}$, ${\rho}_{2}$, and ${\rho}_{\text{d}}$ are the total concentrations of A and B polymers and dimers (measured in polymeric units), that is, ${\rho}_{1}={c}_{1}/{L}_{1}$, ${\rho}_{2}={c}_{2}/{L}_{2}$, and ${\rho}_{\text{d}}={c}_{\text{d1}}/{L}_{1}={c}_{\text{d2}}/{L}_{2}$. Equations (9) and (10) are consistent with the definitions of the dissociation constants of a dimer and of an independent bond, respectively.
The freeenergy density due to nonspecific interactions can in general be written as a power expansion in the concentrations (Semenov and Rubinstein, 1998; De Gennes, 1979),
where the sum is over all the species in the system, including free polymers/stickers, dimers and independent bonds, and ${v}_{ij}$ and ${w}_{ijk}$ are two and threebody interaction parameters. For our simulation system, we derive a specific form of ${F}_{\text{ns}}$ by taking into account that (1) we are interested in the strongbinding regime where the magicratio effect is observed, (2) there is no nonspecific interaction between free polymers of different types in our simulation, and (3) nonspecific interactions are only important at high concentrations. The result is
where v_{b} and w_{b} are the two and threebody interaction parameters for a solution of independent bonds. See Appendix 2 for details of the derivation.
Finally, substituting the conditions Equations (7) and (8) into Equation (6), we obtain the total freeenergy density $F={F}_{\text{ni}}+{F}_{\text{s}}+{F}_{\text{ns}}$,
where ${c}_{\text{d1}}$, ${c}_{\text{d2}}$, and ${c}_{\text{b}}$ are the solutions of Equations (7) and (8). Equations (7), (8) and (13) form a complete set which predicts the freeenergy density of the twocomponent associative polymer system at given total global sticker concentrations, c_{1} and c_{2}, of the two species.
Intuitively, in the strongbinding regime, that is when ${c}_{1},{c}_{2}\gg {K}_{\text{d}},{K}_{\text{b}}$, polymers either associate as dimers or as independent bonds depending on their relative free energies. In the limit that dimers are preferred (${\rho}_{\text{d}}=\mathrm{min}({\rho}_{1},{\rho}_{2})$ and ${c}_{\text{b}}=0$), the contribution from specific interactions is
where $\rho =\mathrm{max}({\rho}_{1},{\rho}_{2})$ is the concentration of the majority species in polymeric units. The terms on the right of Equation (14) reflect, respectively, the freeenergy density due to dimer formation, translational entropy of leftover polymers, and loss of translational entropy of the majority species (in effect, the formation of each dimer removes the translation entropy of one free polymer). In the opposite limit that independent bonds are preferred (${c}_{\text{b}}=\mathrm{min}({c}_{1},{c}_{2})$ and ${\rho}_{\text{d}}=0$),
where $c=\mathrm{max}({c}_{1},{c}_{2})$, and the terms are analogous to those in Equation (14). Numerical studies show that the full ${F}_{\text{s}}({c}_{1},{c}_{2})$ in Equation (6) is always well approximated by the lower of the two limiting values of ${F}_{\text{s}}$ (Equation (14) and (15)).
In which regions of concentration space are dimers versus independent bonds preferred? For a magicnumber system composed of two polymer species of valence L at equal sticker stoichiometry, ${F}_{\text{s}}^{\text{dim}}/{k}_{\text{B}}T=\rho \mathrm{ln}({K}_{\text{d}}e/\rho )$ and ${F}_{\text{s}}^{\text{ind}}/{k}_{\text{B}}T=c\mathrm{ln}({K}_{\text{b}}e/c)$. Comparing the two expressions, dimers are favored at low concentrations, whereas a network of independent bonds is favored at high concentrations. The transition occurs when ${F}_{\text{s}}^{\text{dim}}={F}_{\text{s}}^{\text{ind}}$, that is at concentration ${c}_{0}=e{({K}_{\text{b}}^{L}/({K}_{\text{d}}L))}^{1/(L1)}$. Away from equal stoichiometry, the transition occurs at a lower concentration ${c}_{s}={c}_{0}{(s1)}^{s1}{s}^{s}$, where $s=\mathrm{max}({c}_{1},{c}_{2})/\mathrm{min}({c}_{1},{c}_{2})>1$ (see Appendix 2 for details). As c_{s} decreases rapidly with increasing s (Figure 6A inset, white curve), the preference for dimers over a gel exhibits a sharp peak around equal stoichiometry.
To give a concrete example of the above analysis, we extract the values of ${K}_{\text{d}}$ for dimers from simulations, choose a value of ${K}_{b}$ for independent bonds close to the dissociation constant of a pair of stickers (see Appendix 2 for details), and numerically solve Equation (7) and (8) for ${c}_{\text{d1}}$, ${c}_{\text{d2}}$, and ${c}_{\text{b}}$ to find the fraction of stickers in dimers and independent bonds for all concentrations $({c}_{1},{c}_{2})$. We find that indeed for polymers of equal valence, dimers are favored at low concentrations and independent bonds at high concentrations. The dimer dominated region extends sharply to higher concentrations in a narrow zone around the diagonal, as quantitatively captured by c_{s} (Figure 6A inset and Appendix 2—figure 1A). For polymers of similar but unequal valence, the dimer dominated region extends to higher concentrations along the direction of equal polymer stoichiometry (Figure 6B inset and Appendix 2—figure 1B).
Finally, to extract the binodal phase boundaries, we substitute the values of ${c}_{\text{d1}}$, ${c}_{\text{d2}}$, and ${c}_{\text{b}}$ into Equation (13) to first obtain the free energy as a function of c_{1} and c_{2}. The freeenergy landscape has two basins, one at small concentrations corresponding to the dilute dimerdominated phase, and one at high concentrations corresponding to the dense independentbonddominated gelphase (Appendix 2—figure 2). We locate the phase boundaries by applying convexhull analysis to this freeenergy landscape (see Appendix 2).
Does the dimergel theory capture the magicratio effect revealed by our simulations? Figure 6A and B show the phase diagrams of A_{8}:B_{8} and A_{8}:B_{7} systems. In both cases, the phase boundaries on the dilute side extend sharply into the twophase region along the direction of equal polymer stoichiometry (tie lines along this direction are denoted by black dots). Figure 6C and D show the dilute and densephase concentrations for systems A_{8}:B_{610} at global sticker stoichiometries 8:610. Notably, the dilutephase concentrations are substantially shifted up around the diagonal, verifying the magicratio effect observed in simulations (Figure 5A).
One of the major assumptions of the dimergel theory is a meanfield approximation. Meanfield theory ignores correlations in binding between stickers in the same chain, and therefore has been applied to long chains in the weak binding regime (such that not every sticker is bound) (Prusty et al., 2018; Choi et al., 2020b). Our dimergel theory bypasses this stringent requirement by explicitly assuming the dilutephase components to be dimers, and only considers stickers to associate independently in the dense phase. This approximation captures a key feature of the dense phase, namely that a single polymer binds to multiple partners. Nevertheless, because stickers belonging to the same polymer are tethered together with relatively short linkers in our simulations, correlations in binding exist (Appendix 2—figure 5A). Therefore, what should be considered to be ‘independent’ is not individual stickers but rather segments of the binding correlation length (∼1.8 stickers). The dense phase of a valence 14 system is thus more accurately described by the theory at valence $14/1.8\approx 8$. We therefore present results for valence eight systems in Figure 6. (The theoretical phase diagrams and the dilute and densephase concentrations for valence 14 systems also verify the magicratio effect (Appendix 2—figure 4)).
The dimergel theory has only a handful of parameters: the valences L_{1} and L_{2} of polymers A and B, the dissociation constants ${K}_{\text{d}}$ and ${K}_{\text{b}}$ of dimers and independent bonds, and the nonspecific interaction parameters ${v}_{\text{b}}$ and ${w}_{\text{b}}$. How are the phase boundaries and the magicratio effects determined collectively by these parameters? If valence is increased while keeping all other parameters fixed in the theory, for equal valence polymers we find that the dilutephase concentration decreases, while the densephase concentration increases, and the peak with respect to stoichiometry is enhanced in terms of the dilutephase peaktovalley ratio (Appendix 2—figure 6A and B). If valence is increased for unequal valence polymers, we observe that the shape of the dilute phase boundary transitions from a shoulder to a peak (Appendix 2—figure 6C and D). All these features are consistent with the simulation results in Figure 4.
For the theory to agree quantitatively with the phase boundaries from simulations, we find that smaller values of nonspecific interaction parameters are necessary for higher valence systems (Appendix 2—figure 6C and D). Intuitively, this follows because higher valence polymers have more backbone bonds, which bring bound sticker pairs closer together in the dense phase – effectively reducing the nonspecific repulsion between them. Finally, the dimergel theory also predicts that the magicratio effect disappears in the weakbinding regime (Appendix 2—figure 7), consistent with our simulation results (Figure 3).
Discussion
Intracellular phase separation is driven by multivalent interactions between macromolecules. These interactions are separated into two classes (Ditlev et al., 2018; Pak et al., 2016): (1) specific interactions, such as binding between protein domains, are relatively strong and involve specific partners and (2) nonspecific interactions, such as electrostatic and hydrophobic interactions, which are much weaker, more generic, and nonsaturable. Multivalent systems with specific interactions allow for ‘orthogonal’ condensates to form: the specific interactions holding together one class of droplets will typically not interfere with those holding together another class. Motivated by the key role of specific interactions in intracellular phase separation, we focused on exploring the effects of specific interactions on the phase boundaries of twocomponent associative polymers. Specifically, we combined coarsegrained molecular dynamics simulations and analytical theory to examine the individual roles of valence, stoichiometry, and binding strength on the phase boundaries. In particular, we identified a magicratio effect: for sufficiently strong binding, phase separation is strongly suppressed at equal polymer stoichiometry.
The magicratio effect occurs exclusively in the strongbinding regime. Are specific proteinprotein, proteinRNA, and RNARNA interactions strong enough to lead to the magicratio effect? The onset of the effect in our simulations occurs around ${U}_{0}=9{k}_{\text{B}}T$ (Figure 3A), which corresponds to a stickersticker dissociation constant $K}_{\text{d}}=0.4\phantom{\rule{thinmathspace}{0ex}}\text{mM$. This value is consistent with the onset ${K}_{\text{d}}$ of 1–2.5 mM estimated from 3D lattice simulations with one polymer and one rigid component (Xu et al., 2020). For comparison, the measured ${K}_{\text{d}}$ for a SUMO protein domain with a SIM peptide is 0.01 mM (Banani et al., 2016) and the ${K}_{\text{d}}$ for an SH3 domain and a PRM peptide is 0.35 mM (Li et al., 2012). Thus for systems as strongly interacting as SUMOSIM or SH3PRM, the magicratio effect in principle should manifest in their phase diagrams. However, the magicratio effect has not been observed in these systems (Li et al., 2012; Banani et al., 2016), possibly due to size and linker length mismatch between the two associating polymers. Furthermore, real biological systems are more complex than our simple model. For example, there can be multipletoone binding, multiple components, and the spacers/linkers can also play nontrivial roles (Banjade et al., 2015; Harmon et al., 2017). Currently, the in vivo relevance of the effects explored in this work remains an open question. Magicratio effects could also manifest in other experimental systems, such as nonbiological polymers, DNA origami (Hu and Niemeyer, 2019), or patchy colloid systems (Bianchi et al., 2011). As an inverse problem, the magicratio effect could be exploited to determine the relative valence of associating biomolecules by measuring their phase diagram.
The magicratio effect allows for novel mechanisms of regulation. Chemical modifications, such as phosphorylation or SUMOylation, which change the effective valence of one component into or out of a magicratio condition could shift the phase boundary as a means of condensate regulation. Cells may also have evolved to avoid magic ratios so as to better promote condensate formation. For example, EPYC1 has valence five and Rubisco has valence eight, and the geometry of binding sites on Rubisco and the length of linkers in EPYC1 are such that they disfavor fullybonded RubiscoEPYC1 dimers even at equal polymer stoichiometry, which suppresses the magicratio effect (He et al., 2020). However, active removal of a terminal EPYC1 binding site, for example by phosphorylation (Turkina et al., 2006), would dramatically change the valence ratio to 1:2, which would then favor stable trimer formation, as previously suggested (Freeman Rosenzweig et al., 2017). We hope that our work will stimulate exploration of magicratio effects in both natural and synthetic multivalent, multicomponent systems.
The simulations and theory presented here are aimed at providing conceptual insights into the phase separation of associating polymers that form onetoone specific bonds. Quantitative descriptions of related real systems will likely require additional features, such as details of molecular shape and flexibility, linker lengths, as well as range and type of interactions. For example, while the magicratio effect is robust with respect to the strength of nonspecific interactions and linker length, these variables do strongly influence phase boundaries. The dilutephase concentrations in our simulations are ∼mM, while the reported values for biological systems are typically tens of μM or less. The discrepancy is likely due to different strengths of nonspecific attraction, different length scales of steric replusion between stickers, and/or different lengths and flexibilities of the linkers (Bhandari et al., 2021). Indeed, increasing the nonspecific attraction in our simulations by a small amount $0.07{k}_{\text{B}}T$ leads to a 50% reduction in the dilutephase concentration (Appendix 1—figure 4A). Reducing the steric repulsion between beads of the same type has a similar effect (Appendix 1—figure 5A). More significantly, increasing the mean linker length from 4.7 nm to 5.9 nm leads to a more than 10fold reduction in the dilutephase concentration (Appendix 1—figure 6A). On the other hand, the densephase concentration strongly depends on the steric repulsion — increasing the sticker size from 2.5 to 2.9 nm decreases the dense phase concentration by a factor of 2 (Appendix 1—figure 5B). This is consistent with results from previous studies on the role of linkers: a selfavoiding random coil linker which occupies a large volume can substantially lower the densephase concentration and even prevent phase separation (Harmon et al., 2017). Future work will explore the interplay between specific and nonspecific interactions, and other molecular properties, and their roles in determining the physical properties of droplets, such as surface tension, viscosity, and rate of exchange between phases.
Appendix 1
Modeling twocomponent multivalent associative polymers
We perform coarsegrained moleculardynamics simulations using LAMMPS (Plimpton, 1995) to simulate twocomponent multivalent associative polymers. Individual polymers are modeled as linear chains of spherical particles connected by harmonic bonds (Appendix 1—figure 1A, type A polymer in blue and type B polymer in yellow). Bonds are modeled using a harmonic potential (Appendix 1—figure 1C, left)
where r_{b} = 4.5 nm is the mean bond length, $k=20{k}_{\text{B}}T/{r}_{\text{b}}^{2}$ is the bond stiffness, ${k}_{\text{B}}$ is the Boltzmann constant, and T = 300 k is room temperature.
Stickers of the same type interact through a softened, truncated LennardJones potential (Appendix 1—figure 1C, middle) ^{11}See LAMMPS manual at https://lammps.sandia.gov/doc/Manual.html for details about this potential.
where $\u03f5=0.15{k}_{\text{B}}T$, $\lambda =0.68$, $\sigma =3.5\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$, and $r}_{\text{c}}=5\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m$. These parameters effectively lead to a sticker of diameter $d\simeq 3\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$ and a weak attractive tail of depth $0.06{k}_{\text{B}}T$. The weak attractive tail is employed solely to promote a more compact dense condensate.
Stickers of different types interact through an attractive potential (Appendix 1—figure 1C, right)
where ${U}_{0}=14{k}_{\text{B}}T$ is used in all simulations in the main text, except as indicated in Figure 3. The attraction cutoff distance is $r}_{0}=2\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m$. Note that due to the strong repulsion between stickers of the same type, simultaneous binding of two stickers of one type to a sticker of the other type is energetically highly disfavored (Appendix 1—figure 1D). This ensures onetoone binding of stickers of different types (Appendix 1—figure 1B). Across our simulations, on average the fraction of stickers that have more than one partner is less than 0.001%.
Phase equilibration and data recording
Each system consists of n_{1} and n_{2} polymers of types A and B, respectively. The number of polymers are determined by their valences/lengths (L_{1} and L_{2}) and global A:B sticker stoichiometry (s) through
The round function is used as we can only simulate an integer number of polymers. For all simulations except those in Figure 5 in the main text, we use $N=2500$, so the total number of stickers is around 2500.
Simulations are equilibrated using a Langevin thermostat in the $NVT$ ensemble at T = 300 K in a box of size 250 nm×50 nm×50 nm with periodic boundary conditions, that is the system evolves according to Langevin, 1908:
where ${\overrightarrow{r}}_{i}$ is the coordinate of particle i, m is its mass, $\gamma $ is the friction coefficient, $\overrightarrow{f}$ is random thermal noise, and the energy $U({\overrightarrow{r}}_{1},\mathrm{\dots},{\overrightarrow{r}}_{N})$ contains all interactions between particles, including harmonic bonds, nonspecific, and specific interactions (Equations (16–18)).
To promote phase equilibrium and ensure that only a single dense condensate is formed, we first initialize the simulation by confining polymers in the region $50\phantom{\rule{thinmathspace}{0ex}}\text{nm}<x<50\phantom{\rule{thinmathspace}{0ex}}\text{nm}$. The attractive interaction between A and B stickers (Equations (18)) is gradually switched on from U_{0} = 0 to 14 over 10^{8} time steps. The Langevin thermostat is applied using a damping factor $\tau =m/\gamma =125\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{s}$, step size $dt=2.5\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{s}$, and mass of particle m=3534.3 ag during this time period. These parameters give the particle the right diffusion coefficient $D={k}_{\text{B}}T/(3\pi \eta d)$ for times longer than $\tau $, where $\eta $ is the water viscosity 0.001 kg/m/s and d the diameter of the particle. This annealing procedure leads to the formation of a dense phase close to its equilibrated concentration. The confinement is then removed, and the system is equilibrated for 10^{8} more time steps to allow the formation of dilute phase and relaxation of the dense phase. After these procedures, we switch to smaller τ=10 ns, dt=0.5 ns, and m=282.7 ag for data recording (D remains the same). The system is relaxed for another 10^{8} steps. The relaxation time of the system depends on the stickersticker bond lifetime; to ensure that the dilute and dense phases are in equilibrium, the above choice of relaxation time before recording corresponds to ∼2000 bond lifetimes. We then recorded the positions of all particles every 10^{6} steps for 400 recordings. For each choice of valence and stoichiometry, we performed 10 simulation replicates with different random seeds. We also checked whether there are systematic deviations between the first and second halves of the recorded simulations, and found consistent results between the two halves.
To test the effect of finite size on the phase boundaries, simulations in Figure 5 in the main text are performed with $N=5000$ and a box of size 315 nm×63 nm×63 nm, that is both total number of stickers and box volume are doubled while the total global sticker concentration remains the same (6.64 mM). Procedures for equilibration and data recording are the same (including the initial confinement region) except systems are relaxed for $5\times {10}^{8}$ steps at $dt=0.5\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{s}$ before recording, as the larger system requires a longer relaxation time.
Determining the phase boundaries
To determine the phase boundaries, we need to obtain the concentrations in the dilute and dense phases. The data we recorded for each system contains 4000 snapshots of polymer configurations (10 replicates and 400 time points each). To measure concentrations, we first group polymers into clusters in each snapshot. Connected stickers are grouped into one cluster: two stickers of the same type are connected if they are neighbors in the same polymer, and two stickers of different types are connected if they are within the attraction distance r_{0} = 2 nm. In most of our simulations, in each snapshot, we observe one large cluster which contains most of polymers, and a few to tens of very small clusters (Appendix 1—figure 2). There is a clear gap between the sizes of these large and small clusters. We define the large cluster as the dense phase, and the smaller clusters as constituents of the dilute phase. In cases where the separation between the dense and dilute phases is unclear, we discard the data set. Appendix 1—figure 2 shows the size distribution of clusters pooled over all snapshots. Appendix 1—table 1 lists all the dilutephase components and their total sticker percentages in the A_{14}:B_{14} system at global sticker stoichiometry 1.21.
To find the dilute and densephase concentrations, we calculate the center of mass of the dense cluster for each snapshot, and recenter the simulation box to this center of mass. We then compute the sticker concentration histogram along the x axis with a bin size 1/50 of box length. The resulting concentration profile has high values in the middle corresponding to the densephase concentration, and low values on the two sides corresponding to the dilutephase concentration. The dilute and densephase concentrations are calculated by averaging the concentration profile over the regions ($x\le 100\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$ or $x\ge 100\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$) and ($10\phantom{\rule{thinmathspace}{0ex}}\text{nm}\le x\le 10\phantom{\rule{thinmathspace}{0ex}}\mathrm{n}\mathrm{m}$), respectively.
Testing the effect of finite simulation size on the phase boundaries
To test the effect of finite size on the phase boundaries, we compare the dilute and densephase concentrations of systems with $N=2500$ to systems with $N=5000$ (Appendix 1—table 2). The systems being compared have the same valence, stoichiometry, and total sticker concentration. Simulations with different total numbers of particles show consistent results, suggesting the effect of finite size is minor.
Determining the percolation threshold
Phase transitions in associative polymeric systems can be thought as phase separation aided percolation, that is, the dense phase of an associative polymer system is a percolating network Semenov and Rubinstein, 1998; Harmon et al., 2017; Choi et al., 2020a. What would be the percolation threshold in our associative polymer system if the density remained homogeneous? To answer this question, we determined the percolation threshold of an A_{8}:B_{8} system at a weak binding strength ${U}_{0}=8{k}_{\text{B}}T$, which avoids phase separation. Briefly, for a given sticker concentration $({c}_{\text{A}},{c}_{\text{B}})$, we perform one simulation at ${U}_{0}=8{k}_{\text{B}}T$ and analyze the size of clusters for all snapshots recorded after the system equilibrates. We judge whether polymers are in a sol or gelstate based on the following gelation criterions: First, for each snapshot if the largest cluster contains more than 70% of the stickers and the second largest cluster contains less than 10% of the stickers in the system, we label this snapshot as having a percolating cluster. Second, for a given system if more than 50% of snapshots have a percolating cluster, we label this $({c}_{\text{A}},{c}_{\text{B}})$ point as a ‘gel’ state. The systems do not meet the gelation criterions are labeled as a ‘sol’ state (Appendix 1—figure 3). It is clear that the dilute/dense phases of the A_{8}:B_{8} system formed at a strong binding strength ${U}_{0}=14{k}_{\text{B}}T$ are in a sol/gelstate. By inspection, the same clear dichotomy applies to the other systems we simulated.
Measuring the dissociation constants
The phase boundaries of associative polymers from theory are very sensitive to model parameters. Here, we extract the dimer and sticker dissociation constants ${K}_{\text{d}}$ and ${K}_{\text{b}}$ from simulations. These values are then utilized in the dimergel theory to obtain the freeenergy density landscape and predict the phase boundaries. Our simulations are performed in the strong binding regime. At our chosen binding strength of ${U}_{0}=14{k}_{\text{B}}T$, the long lifetime of each bond means that dimers of valence ≥4 never dissociate in our simulations. This prevents us from directly extracting the dimer dissociation constants for long polymers. We therefore use a reweighting method Frenkel and Smit, 2001 to obtain the dissociation constant for these dimers.
Briefly, we perform simulations with 1 polymer of type A and 1 polymer of type B of valences L_{1} and L_{2} in a cubic box of side 20L_{1} nm with periodic boundary conditions. The systems are equilibrated using a Langevin thermostat in the $NVT$ ensemble, all the parameters are the same as the ones used for recording the data, except here we use ${U}_{0}=7{k}_{\text{B}}T$ which is half of the original value, in order to allow dimers to dissociate.
Theoretically, the dissociation constant of a dimer is defined as:
where $\beta =1/{k}_{\text{B}}T$, and $({\overrightarrow{r}}_{1}^{\text{A}},\mathrm{\dots},{\overrightarrow{r}}_{{L}_{1}}^{\text{A}})$ and $({\overrightarrow{r}}_{1}^{\text{B}},\mathrm{\dots},{\overrightarrow{r}}_{{L}_{2}}^{\text{B}})$ are the coordinates of stickers in polymers A and B. ${U}_{\text{A}}$(${U}_{\text{B}}$) contain all interactions within the polymer A(B), including bond potentials ${U}_{\text{b}}$ (Equations (16)) and nonspecific interactions ${U}_{\text{r}}$ (Equations (17)). ${U}_{\text{AB}}$ contains all the specific interactions between polymers A and B, that is the sum of ${U}_{\text{a}}s$ (Equations (18)). Integration is over the entire volume V. In the numerator, the integration is further confined to the region where ${U}_{\text{AB}}<0$. Note that for a dissociated dimer ${U}_{\text{AB}}=0$.
To link the simulation with the definition of ${K}_{\text{d}}$, we define three variables in the simulation:
where ${E}_{\text{AB}}={U}_{\text{AB}}({\overrightarrow{r}}_{1}^{\text{A}},\mathrm{\dots},{\overrightarrow{r}}_{{L}_{1}}^{\text{A}},{\overrightarrow{r}}_{1}^{\text{B}},\mathrm{\dots},{\overrightarrow{r}}_{{L}_{2}}^{\text{B}})$. We then have,
where $\u27e8{w}_{1}\u27e9$, $\u27e8{w}_{2}\u27e9$, and $\u27e8{w}_{3}\u27e9$ are the mean values of w_{1}, w_{2}, and w_{3} obtained by averaging over a simulation with ${U}_{0}=7{k}_{\text{B}}T$. C and $\stackrel{~}{C}$ are constants and are the same in all three equations. Therefore,
On the other hand, it can be shown that the binding probability $\u27e8{w}_{2}\u27e9$ for ${U}_{0}=7{k}_{\text{B}}T$ is
Combining Equations (27) and (28), we have
Appendix 1—table 3 shows a list of dissociation constants from simulations. The reweighting method provides very accurate stickersticker and dimerdimer dissociation constants, as confirmed by comparing with theory and direct simulations for stickers and polymers of length $L=2$. The dissociation constants obtained by the reweighting method are used in the dimergel theory to predict the phase boundaries (Appendix 2).
Effects of nonspecific interactions and linker length
In the main text, we focused on the effects of binding strength, sticker stoichiometry, and polymer valences on the phase boundaries of twocomponent systems. Here, we explore the effects of nonspecific interactions and linker length. Increasing the nonspecific attraction between stickers of same type leads to decreased/increased dilute/densephase concentrations (Appendix 1—figure 4). The additional attractive interaction is modeled by a cosinesquared potential
The cosinesquared potential is applied together with the softened, truncated LennardJones potential (Equation (17)). Similarly, reducing the range of repulsion between stickers of the same type leads to decreased/increased dilute/densephase concentrations (Appendix 1—figure 5). Here, the repulsive interaction potential between same type of stickers is replaced by the standard repulsive LennardJones potential
The phase boundaries are also sensitive to the linker length. Here, we model the repulsive interactions between same type stickers by Equation (33), the bond between stickers by an expanded FENE potential
and attraction between different types of stickers by Equation (18). Increasing of mean bond length from 4.7 to 5.9 nm leads to a decrease of the dilutephase concentration by more than a factor of 10 (Appendix 1—figure 6).
Appendix 2
Derivation of freeenergy density for nonspecific interactions
The freeenergy density due to nonspecific interactions can be written as a power expansion in the concentrations Semenov and Rubinstein, 1998; De Gennes, 1979,
where the sum is over all the species in the system, including free polymers/stickers, dimers and independent bonds, and ${v}_{ij}$ and ${w}_{ijk}$ are two and threebody interaction parameters.
In the strongbinding regime where the magicratio effect is observed, bound stickers strongly overlap, so the size of a bound pair is almost that of a free sticker. For simplicity, we therefore assume the interactions between dimers to be the same as between free polymers of the same type (denoted as ${v}_{\text{d}}$ and ${w}_{\text{d}}$), and the interactions between independent bonds to be the same as between free stickers of the same type (denoted as ${v}_{\text{b}}$ and ${w}_{\text{b}}$). When independent bonds are preferred (i.e. when ${\rho}_{\text{d}}\approx 0$), the freeenergy density for nonspecific interactions is
where c_{1}, c_{2}, and ${c}_{\text{b}}$ are the total concentrations of polymer A, B, and independent bonds in sticker units. ${c}_{1}{c}_{\text{b}}$ and ${c}_{2}{c}_{\text{b}}$ are therefore the concentrations of free sticker A and B. Note that in our simulations, there is no nonspecific interaction between free polymers of different types. Therefore, all v and w terms involving different free species are 0. Equation (36) simplifies to
Similarly, when dimers are preferred (i.e. when ${c}_{\text{b}}\approx 0$), the freeenergy density for nonspecific interactions is
where ${\rho}_{1}$, ${\rho}_{2}$, and ${\rho}_{\text{d}}$ are the total concentrations of polymer A, B, and dimers in polymeric units.
As nonspecific interactions are only important at high concentrations, we simply set ${F}_{\text{ns}}={F}_{\text{ns}}^{\text{ind}}$. Further, in the strongbinding regime, ${c}_{\text{b}}\approx \mathrm{min}({c}_{1},{c}_{2})$, so
Determining model parameters
The developed dimergel theory has only four parameters for a given system A_{L1}:B_{L2}: the dissociation constants of dimers and independent bonds ${K}_{\text{d}}$ and ${K}_{\text{b}}$, and the nonspecific interaction parameters ${v}_{\text{b}}$ and ${w}_{\text{b}}$. The four parameters together determine the competitiveness of the dilute dimerphase with the dense gelphase.
We have extracted the values of ${K}_{\text{d}}$ from simulations (Appendix 1—table 3). Physically, we expect ${K}_{\text{b}}$ to be close to the stickersticker dissociation constant $5.7\mu \text{M}$ (Appendix 1—table 3). It is not exactly the same because the bonds in the condensate are tethered by the backbones of the polymers. We expect the nonspecific interaction parameters to be approximately ${v}_{\text{b}}=2{B}_{2}$ and ${w}_{\text{b}}=3{B}_{3}$, where B_{2} and B_{3} are the second and third virial coefficients Katsura, 1959: ${v}_{\text{b}}=6.8\times {10}^{2}{\text{mM}}^{1}$ and ${w}_{\text{b}}=2.2\times {10}^{3}{\text{mM}}^{2}$ for hard spheres of diameter 3 nm (the size of a sticker in the simulations), respectively.
The predicted phase boundaries are sensitive to these parameters. We thus tune ${K}_{\text{b}}$, ${v}_{\text{b}}$, and ${w}_{\text{b}}$ around their estimated values to match the dilute and densephase concentrations of the A_{8}:B_{8} system from simulation. Specifically, we yield ${K}_{\text{b}}=3.8\times {10}^{3}\text{mM}$, ${v}_{\text{b}}=9\times {10}^{2}{\text{mM}}^{1}$, and ${w}_{\text{b}}=7\times {10}^{3}{\text{mM}}^{2}$. These parameters are used for A_{8}:B_{610} systems. Further discussions on how model parameters affect phase boundaries can be found in later sections.
Derivation of the transition concentration c_{s}
In the strongbinding regime, the freeenergy density contributions from specific interactions in the dimerdominated and independent bonddominated limit are, respectively,
where ${\rho}_{+}=\mathrm{max}({c}_{1}/{L}_{1},{c}_{2}/{L}_{2})$, ${\rho}_{}=\mathrm{min}({c}_{1}/{L}_{1},{c}_{2}/{L}_{2})$, ${c}_{+}=\mathrm{max}({c}_{1},{c}_{2})$, and ${c}_{}=\mathrm{min}({c}_{1},{c}_{2})$. Equations (40) and (41) are the same as Equations (14) and (15).
At given (c_{1},c_{2}), whether the system will form dimers or independent bonds depends on their relative free energies. For equal valence polymers ${L}_{1}={L}_{2}=L$, letting ${c}_{+}=s{c}_{}$, Equations (40) and (41) become
Comparing the two expressions, dimers are favored at low concentrations (${F}_{\text{s}}^{\text{dim}}<{F}_{\text{s}}^{\text{ind}}$), and independent bonds are favored at high concentrations (${F}_{\text{s}}^{\text{ind}}<{F}_{\text{s}}^{\text{dim}}$). The transition occurs when ${F}_{\text{s}}^{\text{dim}}={F}_{\text{s}}^{\text{ind}}$, i.e. at the concentrations ${c}_{}(s)={c}_{0}{(s1)}^{s1}{s}^{s}$ and correspondingly ${c}_{+}(s)={c}_{0}s{(s1)}^{s1}{s}^{s}$ where ${c}_{0}=e{({K}_{\text{b}}^{L}/({K}_{\text{d}}L))}^{1/(L1)}$. The boundary between dimer and independent bonddominated regions is described by $({c}_{+}(s),{c}_{}(s))$ and $({c}_{}(s),{c}_{+}(s))$, respectively, in the lower and upper halves of the $({c}_{1},{c}_{2})$ plane (Appendix 2–figure 1A, white curve).
Solving reaction Equations (7) and (8)
The high powers in Equation (7) and the small value of ${K}_{\text{d}}$ make it difficult to find numerical solutions of ${c}_{\text{d}}$ and ${c}_{\text{b}}$ accurately. To overcome this difficulty, we define a variable $\lambda ={c}_{1}{c}_{\text{d1}}{c}_{\text{b}}$ when $\rho}_{1}\le {\rho}_{2$ and $\lambda ={c}_{2}{c}_{\text{d2}}{c}_{\text{b}}$ when ${\rho}_{1}>{\rho}_{2}$, and rewrite Equations (7) and (8) in terms of $\lambda $. Specifically,
for $\rho}_{1}\le {\rho}_{2$, and
for ${\rho}_{1}>{\rho}_{2}$. We then solve Equations (44) and (45) using the MatLab function vpasolve with the constraints $0<\lambda <{c}_{1}$ and $0<\lambda <{c}_{2}$, respectively. vpasolve provides all solutions within the specified range. When multiple solutions coexist, we take the one with the lowest ${F}_{\text{s}}$. Numerical solutions of ${\rho}_{\text{d}}$ and ${c}_{\text{b}}$ for A_{8}:B_{8} and A_{8}:B_{7} systems are shown in Appendix 2–figure 1, where the fraction of stickers in dimers is defined as ${\rho}_{\text{d}}/\mathrm{min}({\rho}_{1},{\rho}_{2})$ and fraction of stickers in independent bonds is defined as ${c}_{\text{b}}/\mathrm{min}({c}_{1},{c}_{2})$.
Determining phase boundaries and tie lines
We obtain the freeenergy landscape by substituting the numerical solutions of ${\rho}_{\text{d}}$ and ${c}_{\text{b}}$ into Equation (13) (Appendix 2–figure 2). We locate the phase boundaries by applying convexhull analysis to this freeenergy landscape using the MatLab function convhull (Figure 6A and B).
To find the tie line going through a given initial concentration point $({c}_{1}^{\text{in}},{c}_{2}^{\text{in}})$, we adopt a modified vector method Marcilla Gomis, 2011. We first draw a line though this point along a direction defined by an angle $\alpha $, and then find the crossing points between this line and the phase boundaries, that is $({c}_{1}^{\text{dil}},{c}_{2}^{\text{dil}})$ and $({c}_{1}^{\text{den}},{c}_{2}^{\text{den}})$. The ‘true’ $\alpha $ is the one that minimizes the freeenergy density of mixing $pF({c}_{1}^{\text{dil}},{c}_{2}^{\text{dil}})+(1p)F({c}_{1}^{\text{den}},{c}_{2}^{\text{den}})$, where $p={r}^{\text{den}}/({r}^{\text{dil}}+{r}^{\text{den}})$ is the volume fraction of dilute phase and $1p$ is the volume fraction of dense phase, and ${r}^{\text{dil}}$ and ${r}^{\text{den}}$ are the distances between the point $({c}_{1}^{\text{in}},{c}_{2}^{\text{in}})$ and the two points $({c}_{1}^{\text{dil}},{c}_{2}^{\text{dil}})$ and $({c}_{1}^{\text{den}},{c}_{2}^{\text{den}})$, respectively.
In order to compare with the simulated dilute and densephase concentrations, we use the initial concentrations from simulations for the specified system at given stoichiometry s: ${c}_{1}={c}_{\text{t}}s/(1+s)$ and ${c}_{2}={c}_{\text{t}}/(1+s)$, where the total sticker concentration is c_{t }= 6.64 mM. We then find the tie line going through this initial concentration point and the corresponding $({c}_{1}^{\text{dil}},{c}_{2}^{\text{dil}})$ and $({c}_{1}^{\text{den}},{c}_{2}^{\text{den}})$. For simplicity, we show the total dilute and densephase concentrations ${c}^{\text{dil}}={c}_{1}^{\text{dil}}+{c}_{2}^{\text{dil}}$ and ${c}^{\text{den}}={c}_{1}^{\text{den}}+{c}_{2}^{\text{den}}$ in Figure 6C and D.
A simple approximation for the freeenergy density F
Finding the numerical solutions of Equations (7) and (8) becomes difficult with increasing valence. From the full solutions of ${c}_{\text{d}}$ and ${c}_{\text{b}}$ (Appendix 2–figure 1), we see that in the dimerdominated region ${c}_{\text{b}}$ is almost 0, and in the independent bonddominated region ${c}_{\text{d}}$ is almost 0. Therefore, we can approximate the freeenergy density as the lower value of the two limiting cases
where
and
are solutions of Equations (9) and (10). Equation (46) provides a very good approximation to the full expression for F (Equation (13)). The phase diagrams derived from Equation (13) and (46) are almost identical for A_{14}:B_{14} (Appendix 2–figure 3). Results in Appendix 2–figure 4 are obtained with this approximation (Equations (46–50)) as there are difficulties solving Equations (7) and (8) numerically for A_{14}:B_{1216} systems due to their high valences.
Correlated binding in the dense phase
In the dimergel theory, we assume that stickers of different types can associate independently in the dense phase. However, as stickers belonging to the same polymer are tethered together, neighboring stickers in one polymer are more likely to bind to neighboring stickers in another polymer, that is there are correlations in binding. To quantify this correlation, we first identify consecutive segments in a polymer that bind to consecutive segments in another polymer. (For example, if in polymer 1 of type A, stickers 1, 2, 3, and 4 bind to stickers 2, 4, 3, and 5 of one polymer of type B, sticker 5 binds to sticker 8 of a second polymer of type B, and stickers 6, 7, and 8 bind to stickers 1, 2, and 3 of a third polymer of type B, then there are three individual segments in polymer 1 of type A with lengths 4, 1, and 3.) Clearly, what should be considered to be ‘independent’ are not individual stickers but rather these consecutively bound segments. To quantify the length of these segments, we measure the probability $p(l)$ that a bound sticker is in a segment of length l. Appendix 2—figure 5 shows the probability distribution $p(l)$ for simulated A_{8}:B_{8} and A_{14}:B_{14} systems at equal stoichiometry. The mean length of ‘independent’ segments is 1.8 for both cases.
Effects of model parameters on phase boundaries
The dimergel theory has only a handful of parameters: the valences L_{1} and L_{2} of polymers A and B, the dissociation constants ${K}_{\text{d}}$ and ${K}_{\text{b}}$ of dimers and independent bonds, and the nonspecific interaction parameters ${v}_{\text{b}}$ and ${w}_{\text{b}}$. These parameters together determine the competitiveness of the dilute dimerphase with the dense gelphase. We first fix ${K}_{\text{b}}=3.8\times {10}^{3}\text{mM}$, ${v}_{\text{b}}=9\times {10}^{2}{\text{mM}}^{1}$, and ${w}_{\text{b}}=7\times {10}^{3}{\text{mM}}^{2}$ (see previous section Determining Model Parameters for how these parameters are derived, and note that the values of ${K}_{\text{d}}$ are taken directly from simulations [Appendix 1—table 3]), and explore the dependence of phase boundaries on the valences L_{1}, L_{2} and on the stoichiometry.
Appendix 2–figure 6A and B show phase diagrams and dilutephase concentrations for A_{4}:B_{4} to A_{14}:B_{14} systems. For these equal valence systems, the dilute/densephase concentrations decreases/increases with increasing valence, and the magicratio effect with respect to stoichiometry is enhanced with increasing valence in terms of dilutephase peaktovalley ratio. Appendix 2–figure 6C and D show phase diagrams and dilutephase concentrations for A_{4}:B_{3} to A_{14}:B_{13} systems. Note that the shape of the dilute phase boundary transitions from a shoulder to a peak with increasing valence. All these features are consistent with the simulation results in Figure 4. Appendix 2–figure 7B shows the dilutephase concentrations for A_{8}:B_{610} systems at equal sticker stoichiometry. The dilutephase concentration is sharply peaked at A_{8}:B_{8}, consistent with the simulation results in Figure 3.
However, the dilutephase concentration of A_{14}:B_{14} does not quantitatively agree with simulation results (Appendix 2–figure 6B vs. Figure 4A). Also, the dilutephase concentrations for unequal valence systems do not decrease with increasing valence (Appendix 2–figure 6D vs. Figure 4C). What is the origin of these discrepancies? Intuitively, the dense phase properties are determined by ${K}_{\text{b}}$, ${v}_{\text{b}}$, and, ${w}_{\text{b}}$. Using the same values of these parameters for systems with different valences means that we are treating the dense phases of these systems as exactly equivalent. However, there are more polymer backbone bonds in higher valence systems. These backbone bonds, from a meanfield point of view, act like attractive potentials between stickers, which effectively reduces the nonspecific repulsion between stickers. Therefore, the dense phases of higher valence systems are energetically favored, and we expect correspondingly lower values of ${v}_{\text{b}}$ and ${w}_{\text{b}}$ for valence 14 systems compared to valence 8 systems. Indeed, we find that somewhat smaller nonspecific interaction parameters, ${v}_{\text{b}}=7\times {10}^{2}{\text{mM}}^{1}$ and ${w}_{\text{b}}=5\times {10}^{3}{\text{mM}}^{2}$, lead to quantitative agreement of the dilute and densephase boundaries for A_{14}:B_{14} and A_{14}:B_{13} systems with the simulation results (Appendix 2–figure 6B and D, curves labeled with *, compared to Figure 4A and C). We thus use these parameters for A_{14}:B_{1216} systems in Appendix 2–figure 4.
Finally, the dimergel theory also predicts that the magicratio effect disappears in the weakbinding regime (Appendix 2–figure 7A and B, gray curves), consistent with the simulation results (Figure 3A).
Data availability
Codes for generating data in this manuscript can be found at https://github.com/yaojunz/matlablammpsphasediagram/tree/codes (copy archived at https://archive.softwareheritage.org/swh:1:rev:eff8367b00e1d12c17542bb9d03d85960a3e53e8/).
References

Biomolecular condensates: organizers of cellular biochemistryNature Reviews Molecular Cell Biology 18:285–298.https://doi.org/10.1038/nrm.2017.7

StructureFunction properties in disordered condensatesThe Journal of Physical Chemistry B 125:467–476.https://doi.org/10.1021/acs.jpcb.0c11057

Patchy colloids: state of the art and perspectivesPhysical Chemistry Chemical Physics 13:6397–6410.https://doi.org/10.1039/c0cp02296a

LASSI: a lattice model for simulating phase transitions of multivalent proteinsPLOS Computational Biology 15:e1007028.https://doi.org/10.1371/journal.pcbi.1007028

Physical principles underlying the complex biology of intracellular phase transitionsAnnual Review of Biophysics 49:107–133.https://doi.org/10.1146/annurevbiophys121219081629

Sequence determinants of protein phase behavior from a coarsegrained modelPLOS Computational Biology 14:e1005941.https://doi.org/10.1371/journal.pcbi.1005941

Who's in and who's OutCompositional control of biomolecular condensatesJournal of Molecular Biology 430:4666–4684.https://doi.org/10.1016/j.jmb.2018.08.003

Thermodynamics of high polymer solutionsThe Journal of Chemical Physics 10:51–61.https://doi.org/10.1063/1.1723621

From DNA nanotechnology to material systems engineeringAdvanced Materials 31:1806294.https://doi.org/10.1002/adma.201806294

Solutions of long chain compoundsThe Journal of Chemical Physics 9:440.https://doi.org/10.1063/1.1750930

Fourth virial coefficient for the square well potentialPhysical Review 115:1417–1426.https://doi.org/10.1103/PhysRev.115.1417

Sur la théorie Du mouvement brownien on the theory of brownian motionComptes Rendus De l'Académie Des Sciences 146:530–533.

Fast parallel algorithms for ShortRange molecular dynamicsJournal of Computational Physics 117:1–19.https://doi.org/10.1006/jcph.1995.1039

Thermodynamics of associative polymer blendsMacromolecules 51:5918–5932.https://doi.org/10.1021/acs.macromol.8b00661

Thermoreversible gelation in solutions of associative polymers. 1. staticsMacromolecules 31:1373–1385.https://doi.org/10.1021/ma970616h

The mechanisms of PMLnuclear body formationMolecular Cell 24:331–339.https://doi.org/10.1016/j.molcel.2006.09.013

Rigidity enhances a magicnumber effect in polymer phase separationNature Communications 11:1–8.https://doi.org/10.1038/s41467020153956
Article and author information
Author details
Funding
National Science Foundation (PHY 1734030)
 Yaojun Zhang
 Bin Xu
 Benjamin G Weiner
 Yigal Meir
 Ned S Wingreen
National Science Foundation (PHY 1521553)
 Benjamin G Weiner
 Ned S Wingreen
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Guanhua He, Martin Jonikas, and Daniel Lee for insightful suggestions and comments. This work was supported in part by the National Science Foundation, through the Center for the Physics of Biological Function (PHY1734030), and NSF grant PHY1521553.
Copyright
© 2021, Zhang et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,868
 views

 517
 downloads

 64
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Physics of Living Systems
We propose the Self Returning Excluded Volume (SREV) model for the structure of chromatin based on stochastic rules and physical interactions. The SREV rules of return generate conformationally defined domains observed by singlecell imaging techniques. From nucleosome to chromosome scales, the model captures the overall chromatin organization as a corrugated system, with dense and dilute regions alternating in a manner that resembles the mixing of two disordered bicontinuous phases. This particular organizational topology is a consequence of the multiplicity of interactions and processes occurring in the nuclei, and mimicked by the proposed return rules. Single configuration properties and ensemble averages show a robust agreement between theoretical and experimental results including chromatin volume concentration, contact probability, packing domain identification and size characterization, and packing scaling behavior. Model and experimental results suggest that there is an inherent chromatin organization regardless of the cell character and resistant to an external forcing such as RAD21 degradation.

 Physics of Living Systems
A hallmark of biomolecular condensates formed via liquidliquid phase separation is that they dynamically exchange material with their surroundings, and this process can be crucial to condensate function. Intuitively, the rate of exchange can be limited by the flux from the dilute phase or by the mixing speed in the dense phase. Surprisingly, a recent experiment suggests that exchange can also be limited by the dynamics at the droplet interface, implying the existence of an ‘interface resistance’. Here, we first derive an analytical expression for the timescale of condensate material exchange, which clearly conveys the physical factors controlling exchange dynamics. We then utilize stickerspacer polymer models to show that interface resistance can arise when incident molecules transiently touch the interface without entering the dense phase, i.e., the molecules ‘bounce’ from the interface. Our work provides insight into condensate exchange dynamics, with implications for both natural and synthetic systems.