Rigidity enhances a magic-number effect in polymer phase separation

Cells possess non-membrane-bound bodies, many of which are now understood as phase-separated condensates. One class of such condensates is composed of two polymer species, where each consists of repeated binding sites that interact in a one-to-one fashion with the binding sites of the other polymer. Biologically-motivated modeling revealed that phase separation is suppressed by a “magic-number effect” which occurs if the two polymers can form fully-bonded small oligomers by virtue of the number of binding sites in one polymer being an integer multiple of the number of binding sites of the other. Here we use lattice-model simulations and analytical calculations to show that this magic-number effect can be greatly enhanced if one of the polymer species has a rigid shape that allows for multiple distinct bonding conformations. Moreover, if one species is rigid, the effect is robust over a much greater range of relative concentrations of the two species.


Cluster size in a null-model with no interactions
In the fully-bonded regime, the average cluster sizes in magic-number systems are lower than those in the nonmagic-number systems ( Fig. 1g-l,m,o). However, average cluster sizes in magic-number systems also become very high at high concentrations. This increase of cluster size occurs because at high polymer concentrations the limited free volume for dimers favors condensation. Is the resulting average cluster size larger than one would expect in a noninteracting system due to accidental overlaps? To estimate the effect of accidental overlaps, we consider a null model with no interactions other than the single-occupancy condition for each species. As shown in Fig. 1m,o, the average cluster sizes in the magic-number systems are close to those in this null model, supporting our conclusion that magic-number systems do not experience enhanced clustering due to the specific interaction.

There is no magic-number effect for weak interactions
In the heat maps of the magic-number systems (Fig. 1h,k) there is a "bump" in the phase boundary, i.e., the cluster size initially increases, and then decreases with specific bond energy. This occurs because in the weak-interaction regime thermal fluctuations lead to free binding sites. The presence of these free binding sites precludes the magicnumber effect; indeed, the system resembles non-magic-number cases and, similarly, the phase boundary moves to lower concentration with increasing interaction strength. However, as the interaction strength is increased further, all possible bonds form and the magic-number effect takes over, shifting the phase boundary back to higher concentration. The net result is the observed non-monotonic behavior ("bump") of the phase boundary with increasing specific bond energy for the magic-number systems.

The magnitude of the magic-number effect increases with valence
Valence is a key parameter in both natural and engineered two-component systems that form specific bonds [1][2][3][4] . It is therefore natural to ask how valence influences the magic-number effect. To theoretically address this question, we studied systems composed of two species of flexible polymers A n :B m over a range of valences ( Supplementary Fig. 2). We observed that non-magic-number systems (A n−1 :B n ) display more clustering than the corresponding magic-number systems A n :B n , despite the higher valence of the latter. Specifically, average cluster size at low concentrations decreases by a factor of 2 from A 3 :B 4 to A 4 :B 4 , 5 from A 5 :B 6 to A 6 :B 6 , and 10 from A 7 :B 8 to A 8 :B 8 ( Supplementary Fig. 2 insets). In Supplementary Fig. 2, we also compared the phase boundaries, i.e., the concentrations at which large clusters begin to form (details in following paragraphs). Note that for each pair, the valence is lower for the nonmagic-number system, which competes with the magic-number effect. Indeed, for the lowest-valence case, A 3 :B 4 has a similar phase boundary to A 4 :B 4 ( Supplementary Fig. 2a). By contrast, for the higher-valence systems, the phase boundaries in the regime of strong interactions occur at substantially higher concentrations for the magic-number systems ( Supplementary Fig. 2b,c). This confirms that higher valence increases the magnitude of the magic-number effect. Note that at lower interaction energies, the presence of unbonded monomers means magic-and non-magicnumber systems behave similarly, which accounts for the "bump" in the phase boundary for A 6 :B 6 and A 8 :B 8 .
We determined approximate phase diagrams for the A 3 :B 4 , A 4 :B 4 , A 5 :B 6 , A 6 :B 6 , A 7 :B 8 , and A 8 :B 8 systems based on fluctuations of cluster size. For each system and concentration, we first increased the non-specific bond from 0 to 0.1 k B T in 0.005 k B T increments, keeping the specific bond energy at 0 k B T . Then the specific bond energy was increased from 0 to 11 k B T in 0.04 k B T increments, while the non-specific bond energy was kept at 0.1 k B T . Each step of annealing was simulated with at least 50,000 Monte-Carlo steps (i.e. proposed moves) per monomer. While increasing the specific bond energy, we recorded the average cluster sizes when the specific bond energy was an integral multiple of 0.4 k B T . This simulation was independently repeated 100 times to compute cluster-size fluctuations, i.e. the standard deviation over 100 simulations of the average cluster size in each simulation.
For each specific bond energy, we expect the peak of cluster-size fluctuations to closely indicate the phase boundary. Therefore, we fit the standard deviation of average cluster size versus concentration with a Gaussian and took the peak to be the concentration corresponding to the phase transition. Cluster-size fluctuations corresponding to different specific energies were fitted independently. The smoothed curves in Supplementary Fig. 2 were obtained via cubic spline, with additional fictitious data points included in the fit (but not shown) to induce a vertical slope for the highest specific bond energies, consistent with the saturation of all bonds.

Three-dimensional model
Simulations were performed on a cubic lattice of 20 × 20 × 20 grid points (or "sites") with periodic boundary conditions. In the model, polymers with different shapes occupy several connected (nearest neighbor) sites such that each monomer occupies one site. There are two species of polymer in each simulation, denoted as A n or B n , with n being the number of monomers in one polymer. A monomer of A and a monomer of B form a specific bond when they occupy the same site in the 3D lattice; no more than one A-monomer and no more than one B-monomer can occupy a site.
The A polymers are considered to be "flexible," such that any configuration of connected nearest-neighbor sites is allowed. The B polymers (B 8R ) are rigid 2 × 2 × 2 cubes.
Systems with only these specific interactions have a very weak effective surface tension between condensed and dilute phases, which prevents formation of dense droplets. Motivated by the existence of weak non-specific interactions between polymers (e.g., due to hydrophobicity), we add a small non-specific interaction between all nearest neighbor monomers as described below, which increases the surface tension between phases and results in denser droplets.
We performed Markov-Chain Monte-Carlo simulations using the Metropolis algorithm. Briefly, in each simulation step we randomly propose a move of the configuration. The move is always accepted if it reduces system energy, and accepted with probability e −(E f −Ei)/kBT , where E f and E i are the final and initial energies, if the move increases system energy. Three categories of moves are proposed: single-flexible-polymer moves, single-rigid moves, and twospecies joint moves. Single-flexible-polymer moves are standard lattice-polymer local moves: the end-point move, the corner move, and the reptation move. Single-rigid moves consist of one-step translations in the six cardinal directions. In the regime of strong specific bonds, the two species are typically held together by multiple specific bonds, which leads to dynamical freezing. To enable the system to better explore configuration space, we include a two-species joint move such that connected clusters of polymers are translated together, without breaking any specific bonds. To obtain thermalized ensembles, we follow a simulated annealing procedure: we keep k B T constant and gradually increase bond strength. We increase the specific bond energy in 0.1 k B T increments, from 0 to the final bond strength (2, 3, 4, or 5 k B T ), while the non-specific bond energy is kept at 0.1 k B T . The system thermalizes over 15,000 Monte-Carlo steps (i.e. proposed moves) per monomer, then results are averaged over the subsequent 15,000 steps. Supplementary Fig. 3 demonstrates that the magic-number effect is present in three dimensions. Snapshots ( Supplementary Fig. 3a-c) show the magic-number system A 8 :B 8R dominated by dimers, whereas the non-magicnumber cases A 7 :B 8R and A 9 :B 8R form large clusters at the same concentration and binding energy. Heatmaps of mean cluster size confirm that the formation of large clusters is suppressed in the magic-number system (Supplementary

Phase separation versus percolation
The biological condensates under consideration are phase-separated droplets and behave as liquids. However, some systems where polymers bond and form connected clusters undergo a homogeneous sol-gel transition. Gels have different physical properties than liquids, which would presumably have implications for biological function. How can we distinguish between phase separation and homogeneous gelation in our two-component multivalent systems?
Within our mean-field theory, the system clearly undergoes phase separation as a first-order phase transition: the free energy is non-convex and the order parameter (density) has a discontinuity, both hallmarks of a first-order transition. Specifically, the non-convex profile of the free energy implies that the ground state of the system in the thermodynamic limit necessarily consists of two phases for all concentrations lying between the points where the tie line is tangent to the free-energy profile, and the densities of the two phases are those corresponding the points of tangency (cf. Supplementary Fig. 8). (Note that although we use mean cluster size as an order parameter throughout this work, here we use density because it exhibits more spatial variation in the finite lattice simulations and facilitates theoretical calculations.) In the Monte Carlo simulations, local density analysis suggests that the macroscopic clusters are phaseseparated droplets rather than percolation clusters arising from homogeneous gelation. First, visual inspection of the interacting system ( Supplementary Fig. 4a) reveals that as the system concentration increases, the cluster or clusters become larger but not denser. In the noninteracting system, in which "bonds" between polymers of different types arise purely from accidental overlaps, this is not the case (Supplementary Fig. 4b). Instead, these systems form large clusters as a simple consequence of percolation. Supplementary Fig. 4c quantifies this difference in terms of the average local density inside the cluster at concentrations where macroscopic clusters are present. To define a local density, we calculate the number of neighbors for every monomer inside a cluster that includes more than half the proteins. (Note that because every site can be occupied by two monomers, the density ranges from 0 to 2.) This quantity is averaged over monomers, and then over Monte Carlo samples. In the interacting system, the local density is quite stable as the system concentration grows. By contrast, the noninteracting ("gel-like") percolation clusters become progressively denser as the concentration grows. This confirms the visual impression from snapshots that the interacting system undergoes phase separation rather than homogeneous gelation. As the macroscopic clusters we observe are always internally connected by a network of specific bonds, we conclude that our system represents a case of gelation driven by phase separation 5,6 .
Finally, it should be noted that suddenly increasing the bond strength (a"quench") can trap the system in a kinetically-arrested state far from equilibrium. This metastable state resembles a network of fibers ( Supplementary  Fig. 5a), which is easily distinguished from the equilibrium droplets ( Supplementary Fig. 5b) by visual inspection.
Supplementary Figure 4. Clusters are phase-separated droplets rather than artifacts of percolation. a Snapshots of simulations of A7:B8R where cluster formation is driven by specific bonds. Parameters: specific bond energy = 10 kBT , non-specific bond energy = 0.1 kBT , A:B monomer ratio = 1. b Snapshots of simulations of A7:B8R where cluster formation is due to random percolation. Parameters: specific bond energy = 0 kBT , non-specific bond energy = 0 kBT , A:B monomer ratio = 1. c The local density inside macroscopic clusters, averaged over all monomers in the cluster then averaged over independent Monte Carlo samples (mean ± SD of Monte Carlo samples).
Clusters are phase-separated droplets at equilibrium rather than kinetically-arrested states. a Snapshot of simulation of A7:B8R after a temperature quench where specific bond energy increased from 0 → 5 kBT . The flexible polymers (An) are red, the rigid cubes (B8R) are blue, and specific bonds where these overlap are green. Parameters: monomer concentration of each species=0.15, non-specific bond energy = 0.1 kBT , A:B monomer ratio = 1. b Snapshot of simulation with same parameters as a, but after simulated annealing to a specific bond energy of 5 kBT .

The role of non-specific interactions
The magic-number effect is due to a competition between translational and conformational entropy in the regime of strong specific bond formation. However, our simulations also include non-specific interactions, which lead to more realistic droplets with a higher density and surface tension (see Methods). How do such non-specific interactions influence the magic-number effect? First, Supplementary Fig. 6 shows that non-specific interactions alone do not lead to a magic-number effect. With the energy of specific bonds set to zero, no large clusters form, and the magic number case A 8 :B 8R (Supplementary Fig. 6b,e) shows no suppression in cluster size compared to the non-magic number cases A 7 :B 8R and A 9 :B 8R . Supplementary Fig. 7 shows the dependence of cluster formation on non-specific interactions for a range of specific bond energies. Weakening the non-specific interactions has two effects: 1) the system must reach higher concentrations before large clusters form and 2) the clusters are less dense. However, this does not change the tendency of magic-number systems to form dimers rather than large clusters. Indeed, the mean cluster size as a function of concentration reveals a magic-number effect even in the absence of non-specific interactions ( Supplementary Fig. 7c). Taken together, Supplementary Fig. 6 and 7 demonstrate that non-specific interactions are neither sufficient nor necessary for the existence of the magic-number effect.

Flory-Huggins theory: a brief review
Flory-Huggins theory 7 provides a simple analytical treatment of phase separation in a polymer-solvent system. The theory estimates both the configurational entropy of the polymers and the enthalpy of polymer-solvent interaction within a mean-field approximation. For strong enough polymer-solvent repulsion, the resulting free-energy density is a non-convex function of the concentration, which implies phase separation.
To briefly review, the Flory-Huggins model considers a lattice on which solvents and polymers occupy sites. Each solvent molecule or monomer occupies one site and all sites are taken to be singly occupied. We consider a system with N p polymers, each with L monomers, and N s solvent molecules. On the lattice with M sites in total, each site has z neighbors.
In order to compute the configurational entropy, we need to estimate the number of configurations in the fully mixed state. Assuming all solvent molecules are identical, they make no contribution to the number of configurations. (Assuming solvent molecules to be distinguishable only contributes a factorial constant, which does not affect the final phase diagram.) Amongst the M lattice sites, we first select N p points on the lattice to be the "head" of each polymer, thus the number of head configurations is: For each of the other monomers of the polymer ("body"), it can choose among (z − 1) sites if the sites are unoccupied. Given the single-occupancy constraint, the number of positions that a "body" monomer can be chosen to occupy is estimated as (z − 1) M −Nocc M , where N occ is the total number of occupied lattice sites before this monomer is placed. Then the number of configurations for all "body" monomers is: Thus the total number of configurations is: For simplicity, it is conventional in Flory-Huggins theory to express the free energy of the well-mixed state relative to a state in which the polymers and solvents are fully separated. In this reference state, the number of configurations is given by Eq. (S3) with M = N p L: The fractional increase of the number of configurations with respect to the reference state is Using Stirling's approximation, we find the entropy change due to mixing: where c = N p L/M is the monomer concentration. Within mean-field theory, the change of enthalpy density, or equivalently internal energy density, upon mixing is: where χ is an interaction parameter, and we express energies in units of the thermal energy k B T . The change in free-energy density is thus This Flory-Huggins free-energy density is a non-convex function of concentration c if χ is large. The total free energy in the non-convex region is minimized by the system separating into two coexisting phases with concentrations at the two common tangent points of a line touching the free-energy curve from below.

Dimer versus condensate theory for magic-number systems
Here, we generalize the Flory-Huggins theory to systems with two polymer species in addition to the solvent. We first focus on the magic-number condition in which the polymers have the same valency L, and restrict our attention to the case with equal monomer concentrations, with each species having N p polymers. In the regime of strong specific bonds between the two polymers, all monomers of the two species are considered to be in bonds, and the internal energy is simply a fixed constant, which we neglect. We therefore consider each site in the lattice to be occupied either by two monomers, one from each polymer species, or by solvent.
i. System with two flexible polymer species (A L :B L ) We model the two-species system separately in two regimes: the high concentration ("dense") regime, dominated by a condensate, and the low concentration ("dilute") regime, dominated by dimers.
In the dense regime, the number of configurations of the fully-bonded system is the product of the number of configurations of the two independent species, times the probability that the subset of sites occupied by the monomers of one species exactly matches the subset of sites occupied by monomers of the other species: where W 1 and W 2 are the total number of lattice configurations for each of the two species considered separately. To compute W 1 and W 2 , we slightly modify Eq. (S3): the translational degrees of freedom are still assumed to contribute while the number of different self-avoiding polymer conformations with a given head position, aka the polymer conformation factor C L , is exactly computed numerically. The total number of configurations is therefore Numerical values of some C L s are: We adopt a mean-field theory for the matching probability, i.e., we consider one species to occupy its complete subset of sites, and then take the probability for the i-th monomer of the other species to fall within the same subset of sites to be P (i) match =

Sites in the matching region & not already occupied
Sites not already occupied = so that, The number of configurations in the dense case is therefore: (S14) We find the free-energy density relative to the state of zero entropy to be: In the dilute regime, we model polymers as fully bonded dimers. The translational degrees of freedom can be approximated by the same expression as in Eq. (S10) now applied to dimers as effectively a single species. Regarding polymer conformations within a dimer, we can exactly numerically compute the number of different ways of forming one dimer with a given head position of one of its polymer types, aka the dimer degeneracy D L . (Note that D L = j d 2 L,j , where d L,j is the number of conformations of a single polymer of length L occupying a specific set j of lattice sites.) The total number of configurations is and the free-energy density relative to a state of zero entropy is Numerical values of some D L s are: To obtain the free-energy density at an arbitrary concentration c, we simply take the minimum of the dense and dilute estimates of the free energy, i.e., (S18) ii. Systems with one flexible polymer species and one rigid species A 8 :B 8R and A 8 :B 8U When one species is a rigid shape, the dimer-versus-condensate theory still applies, but with several modifications. One modification is that the polymer that forms the rigid shape has limited flexibility. Therefore, for W 2 in Eq. (S11) in the dense regime, we replace C 8 = 2172 by C 8R = 4, yielding In the dilute regime for the A 8 :B 8R/U system, the number of ways to form dimers with a given center of mass is reduced from D 8 = 9960 to D 8R = 112 and D 8U = 8 in Eq. (S16), so that Oligomer-versus-condensate theory for non-magic-number systems (A L1 :B L2 ) For systems with two flexible species of polymers with different valencies L 1 and L 2 , we can similarly extend the dimer-versus-condensate theory. In the dilute regime, as the system can no longer form fully saturated dimers, it forms fully bonded oligomers composed of multiple polymers of each species. For simplicity, we consider the case when L 1 and L 2 are co-prime, which includes our simulations with L 1 = L 2 − 1. The smallest oligomer in this case has L 1 L 2 monomers of each species. We focus on systems with equal monomer concentrations in the strongly interacting limit, so that all monomers form specific bonds.
The dense regime is similar to that in the magic-number system, but with the two species having different valencies L 1 and L 2 . The number of configurations is where, similar to Eq. (S11), Following Eq. (S13), we have where c = Np1L1 M = Np2L2 M is the concentration of each monomer type. Then the free-energy density is: In the dilute regime, the system is dominated by oligomers that occupy L 1 L 2 lattice sites. Since an oligomer is typically much larger than a dimer in the magic-number systems, it is impractical to exactly numerically compute the oligomer degeneracy, i.e., the number of ways to form a fully-bonded oligomer with N 1 = L 2 polymers of the first species and N 2 = L 1 polymers of the second species. As a rough approximation, we assume that the dominating configurations form the simplest compact shape, i.e., an L 1 × L 2 rectangle. We then use Eq. (S11) to estimate the number of configurations of each species within an L 1 × L 2 rectangle (so that M → L 1 L 2 ): The resulting expression for the total number of configurations in the dilute case is then similar to Eq. (S16).
and the free energy is