Magic numbers in polymer phase separation – the importance of being rigid

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. Previous modeling of one such two-component system – the pyrenoid, an organelle found in algae and plants – surprisingly revealed that phase separation is suppressed for certain combinations of numbers of binding sites. This phenomenon, dubbed the “magic-number” effect, 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. Our findings advance our understanding of the fundamental physics of two-component polymer-based phase-separation and suggest implications for biological and synthetic systems.


Introduction
In addition to membrane-bound organelles, cells possess non-membrane-bound bodies including nucleoli, Pbodies, and stress granules, which are now understood as phase-separated condensates [1][2][3][4][5][6][7] . Typically, the components of these condensates have a high rate of exchange with the surrounding medium and the condensates themselves are dynamic, rapidly assembling and disassembling in response to specific stimuli [8][9][10] . However, the physical mechanisms responsible for the phase separation of biomolecular condensates remain poorly understood.
Some condensates that offer a tractable model for mechanistic understanding are those composed of two species of multivalent polymers with specific interactions that drive phase separation [11][12][13] . In the simplest case, each polymer consists of repeated domains that interact in a one-to-one fashion with the domains of the other polymer. Such two-component multivalent condensates have been observed in several natural and engineered contexts. One example, the algal pyrenoid, is a carbonfixation organelle, in which the two components essential for assembly are the rigid oligomer Rubisco (the active enzyme in CO 2 fixation) and EPYC1, an unstructured linker protein 14 . Another multivalent system, PML (promyelocytic leukemia) nuclear bodies that repair DNA damage 15 , relies on the Small Ubiquitin-like Modifier (SUMO) domain that interacts with the SUMO Interacting Motif (SIM) to form droplets [16][17][18] . This system inspired in vitro experiments with engineered poly-SUMO and polySIM of various valencies 13 , and, indeed, phase separation was observed with increasing concentrations of the two polymers. Other in vitro two-component systems, e.g., an engineered polySH3-Proline-Rich Motif system 11 and a PTB-RNA system 12 , also form droplets with increasing concentrations.
A striking theoretical prediction regarding these twocomponent multivalent systems is that in the regime of strong binding they will be extremely sensitive to the relative valency of the two species of polymers. Higher valency is normally expected to boost condensate formation 11 . In two-component systems, however, an exception occurs when the valency of one polymer species equals or is an integral multiple of the valency of the other species. In this case, condensation is strongly suppressed in favor of small fully bonded oligomers -a "magic-number" effect reminiscent of the exact filling of atomic shells leading to the unreactive noble gases. For example, in the aforementioned pyrenoid, EPYC1s have 4 nearly identical repeats and Rubisco has 8-fold symmetry. Simulations thus predicted that the magic-number condition 4 + 4 = 8 would lead to preferential formation of trimers composed of two flexible EPYC1s and one rigid Rubisco 10 . Here, we demonstrate that the magic number effect still occurs if all polymers are flexible, but the effect is maximized if one of the polymers is rigid and compact, which entropically favors fully-bonded oligomers over the condensate. Strikingly, while the magic-number effect requires rather precise global stoichiometry when all polymers are flexible, if one of the polymers is rigid the effect occurs over a broad range of stoichiometries, implying that modulation of binding valency can robustly control condensation both in and out of cells.

Results
The magic-number effect with one rigid and one flexible component To examine the role of rigidity in the magic-number effect, we begin with a flexible polymer, motivated by EPYC1, together with a rigid shape, motivated by Rubisco. The flexible polymer with n binding sites is denoted as A n and the rigid 4 × 2 shape is denoted as B 8R ; each binding site on each species is considered to be a "monomer". Previous results which established the existence of a magic-number effect can then be summarized as follows: the A 4 :B 8R system forms stable trimers and thus, compared to the A 3 :B 8R or A 5 :B 8R systems, requires a substantially higher total monomer concentration for the formation of large clusters 10 .
Based on these observations, we predicted that flexible polymers of length 8 together with rigid 4 × 2 shapes should also constitute a magic-number system, in this case based on stable heterodimers, since each constituent has 8 binding sites. To test this idea, we simulated A 7 :B 8R , A 8 :B 8R , and A 9 :B 8R systems (Fig. 1a-c). The heat maps of the average cluster size for different concentrations and specific bond strengths ( Fig. 1g-i) define approximate phase diagrams: regions with smaller cluster sizes correspond to a homogeneous dissolved phase, while regions with larger cluster sizes correspond to an inhomogeneous condensed phase, i.e., a phase-separated regime. As anticipated, the phase boundary for A 8 :B 8R occurs at substantially higher concentration than for A 7 :B 8R or A 9 :B 8R , reflecting the stability of dimers composed of one flexible A 8 and one rigid B 8R polymer ( Fig. 1g-i,m). Note that we define a "dimer" as consisting of one polymer of each species which only form specific bonds with each other.
In Fig. 1g-i, as specific interactions increase from zero, cluster sizes initially increase, but stop increasing above ∼ 7k B T . This behavior is due to saturation of the specific bonds, i.e., for large enough specific bond energies essentially all possible bonds are formed. In this fully bonded regime, as polymer concentration increases, the system undergoes a transition from a uniform dissolved phase to a phase-separated non-uniform "droplet" phase.
What drives phase separation in these two-component systems? The transition is not driven by a competition between entropy and energy as in Flory-Huggins theory 19 , but rather by a competition between different types of entropy. For example, in the droplet phase of the A 8 :B 8R system, free dimers coexist with large gellike clusters. Each dimer has high translational entropy, while each polymer in a cluster has very limited translational entropy. However, this difference in translational entropies is offset by an opposing difference in conformational entropies: in each dimer, the polymer of one species must match all the binding sites of the polymer of the other species, leading to a reduced total conformational entropy. By contrast, the polymers in a large condensed cluster are more independent, binding to multiple polymers of the other species and enjoying a relatively high conformational entropy. Because the importance of translational entropy relative to conformational entropy is reduced as concentration increases, the system transi-tions from a uniform dimer phase to a droplet phase with increasing polymer concentration.
By contrast, in non-magic-number systems (e.g., A 7 :B 8R or A 9 :B 8R ), the polymers cannot form free dimers without unsatisfied binding sites, and instead tend to form large fully-bonded clusters even at low concentrations. This key difference leads to the drastic difference in clustering between magic-number and non-magic-number systems.
To verify that in the A 8 :B 8R magic-number system the phase transition in the fully-bonded regime is driven by a competition between dimers and a condensed phase, we measured the fraction of polymers in the dimer form (Fig. 1n). For magic-number systems in the low concentration regime, the polymers are essentially all in dimers, while at very high concentration the fraction of dimers is low since the system is dominated by a condensate. Notably, the fraction of dimers is always higher than in a null model (see Supplementary Note for details). This confirms the picture of an entropy-driven phase transition in which the A 8 :B 8R magic-number system is dominated by dimers at low concentration and a condensed phase at high concentration. By contrast, the A 7 :B 8R and A 9 :B 8R non-magic-number systems never have a high dimer fraction (Fig. 1n). Even at low concentrations, the polymers tend to aggregate into clusters to form as many bonds as possible.
The magic-number effect is attenuated if all polymers are flexible In the above, we considered one flexible component and one rigid component. In principle, the magic-number effect should still occur if both components are flexible. So does the rigidity of one component matter? To answer this question, we simulated systems of two flexible polymer species, where one species has 7, 8, or 9 binding sites and the other has 8 binding sites (systems denoted as A 7 :B 8 , A 8 :B 8 , and A 9 :B 8 ). Snapshots reveal fewer clusters in A 8 :B 8 than in A 7 :B 8 or A 9 :B 8 ( Fig. 1d-f ), and heat maps of average cluster size confirm that higher concentrations are required for clustering in the A 8 :B 8 system ( Fig. 1j-l). Correspondingly, the A 8 :B 8 system has smaller average cluster sizes (Fig. 1o) and a higher fraction of dimers (Fig. 1p) than the non-magic number systems in the strong-interaction limit. However, compared to the case with one rigid component the magicnumber effect is noticeably weaker when both polymers are flexible: the onset of phase separation occurs at a lower concentration and average cluster sizes are in general larger. What is the origin of this difference?
The conformational entropy of polymer dimers strongly influences clustering A possible reason for the enhanced magic-number effect in systems with one rigid component lies in the number of distinct ways of forming dimers. Given the 4 × 2 rectangular shape of B 8R , an A 8 polymer has 28 different ways to pack inside the rectangle, i.e., the dimer degen- eracy is 4 × 28 = 112, where the factor of 4 comes from the four possible square-lattice orientations of a rectangle with a defined "head" site. By comparison, the dimer degeneracy for A 8 :B 8 is higher, since there are 9960 distinct ways for two length L = 8 polymers to occupy the same lattice sites, with the head of one of the polymers at a defined site. However, to gauge the effect of dimer degeneracy on clustering, one should actually compare the dimer degeneracy to the number of possible conformations of the two polymers considered separately, since the two polymers can sample conformations independently within a condensed cluster. For the A 8 :B 8R system, the number of independent conformations is 4 × 2172, since 2172 is the number of conformations of a polymer of length L = 8 with a defined head site; correspondingly, for the A 8 :B 8 system the number of independent conformations is 2172 × 2172. Therefore, the fold-reduction in the number of conformations upon dimerization in the semi-rigid A 8 :B 8R system is (4 × 2172)/112 = 78, while the fold-reduction for the flexible A 8 :B 8 system is sub-stantially larger (2172 × 2172)/9960 = 474. The greater loss of conformations upon dimerization in the fully flexible A 8 :B 8 system will favor clusters over dimers, and could therefore explain the larger average cluster size for To test this idea, we designed a system with a minimal dimer degeneracy, but otherwise as similar as possible to the semi-rigid A 8 :B 8R system. Namely, we considered a rigid 4 × 2 shape together with a flexible polymer, but with a rule that allows only a single "U-shaped" conformation of the flexible polymer to achieve the full dimer binding energy, while still permitting high conformational entropy of the flexible polymers within large clusters (see Methods). Dimers in this A 8 :B 8U system have an effective degeneracy of only 4 × 2 = 8. As shown in Fig. 2a,b, snapshots of simulations of the A 8 :B 8U system with monomer concentration 0.4 (i.e. 40% of the lattice sites occupied by each species) dramatically confirm that lower dimer degeneracy leads to more clustering. The A 8 :B 8R system with dimer degeneracy 112 is Lower conformational entropy of dimers favors clustering. Snapshots of simulations of a A8:B8R, and b A8:B8U systems. B8U denotes rigid 4 × 2 rectangles that restrict a dimerizing polymer partner to adopt a specific U-shape (see Methods). Parameters: specific bond energy = 10 kBT , non-specific bond energy = 0.1 kBT , A:B monomer ratio = 1, monomer concentration = 0.4. Inset in b shows a fully-bonded U-shaped dimer. c Average cluster size versus monomer concentration for systems in a,b.
mostly composed of dimers, with only a few small clusters, while the A 8 :B 8U system with dimer degeneracy 8 is dominated by one large cluster. Averaged results obtained at different concentrations ( Fig. 2c) confirm that lower dimer degeneracy leads to more clustering.

Relative concentration of monomers also strongly influences clustering
So far, we have considered cases where the two polymer species have the same total number of monomers in the simulation domain, i.e., the same monomer concentration. How do differences in monomer concentration influence clustering for magic-number and non-magicnumber systems? To address this question, we compared the systems A 7 :B 8 to A 8 :B 8 (both polymers flexible) and A 7 :B 8R to A 8 :B 8R (one flexible polymer and one rigid 4 × 2 shape) over a range of relative monomer concentrations. To avoid the confounding effect of changing the total monomer concentration, we fixed the sum of the concentrations of monomers of the two species, while varying the ratio of the two concentrations from 1:2 to 2:1 (Fig. 4).
Both non-magic-number cases, A 7 :B 8 and A 7 :B 8R , display a peak of average cluster size around equal monomer concentration. This is because at equal concentration there are no excess monomers of either species, so the system is driven to form large clusters in order to maximize the total number of specific bonds.
By contrast, the magic-number system A 8 :B 8 with both polymers flexible displays a sharp dip in average cluster size around equal monomer concentration. For strong specific interactions, equal concentrations implies that all monomers are in specific bonds. In this case, for the concentration of 0.3 shown in Fig. 4a, entropic considerations favor dimers over clusters. However, away from equal concentrations the resulting excess of monomers of one polymer species greatly increases the internal conformational entropy of clusters, leading to a rapid increase of average cluster size. The average cluster size is symmetric around equal monomer concentration because the A 8 and B 8 polymers are equivalent.
Strikingly, as shown in Fig. 4b, the average cluster size in the A 8 :B 8R system with one flexible polymer and one rigid shape is asymmetric around equal concentration. Specifically, the system with an excess of the rigid shape B 8R (left of center) is much more clustered than the system with an excess of the flexible polymer A 8 (right of center). An intuitive argument explains this asymmetry: adding an excess rigid shape to a cluster allows many new configurations of the flexible polymers, while adding an excess flexible polymer to a cluster does not allow many new configurations of the rigid shapes. Importantly, this asymmetry implies a robustness of the magic-number effect for systems in which flexible polymers interact with more compact multivalent objects; specifically, a magicnumber ratio of binding sites disfavors clustering -with no fine tuning of concentration -provided the flexible polymer is in excess.

Mean-field theory for the strong interaction limit
Can we go beyond simulations to better understand the role of rigidity in the magic-number effect? The well-known Flory-Huggins theory 19 models the competition between mixing entropy, which favors a dissolved phase, and interactions, which typically favor phase separation. While successful in explaining phase separation in polymer solutions, Flory-Huggins theory only considers non-specific interactions. More recently, Seminov & Rubinstein 20 presented a sticker-solution theory that incorporates specific one-to-one interactions among "stickers" on a single polymer species. The one-component sticker theory predicts both a gelation transition and condensate formation. However, two-component multivalent systems with specific interactions present additional mechanisms and a wider range of behaviors. Notably, as described above, in the limit of high binding energies in these systems the internal energy is constant -essentially all specific bonds are formed -so the phase transition is a purely entropic effect driven by a competition between a dimer (or small oligomer) phase with higher translational entropy and a condensate phase with higher conforma- tional entropy 10 .
To analyze the role of rigidity in the regime of high binding energies, we developed a simple theoretical model for our lattice-based polymer system. Specifically, we consider two species of lattice polymers with equal monomer concentrations in the strong interaction limit. Thus every monomer is in a specific bond.
For one-component lattice-polymer systems, Flory-Huggins theory provides a theoretical framework to compute the configurational entropy and internal energy, and thus the free energy 19 . For a two-component system in the fully-bonded regime, the total internal energy is a constant, and so we focus on the configurational entropy. We approximate the configurational entropy in two limits: the dilute limit dominated by dimers or other oligomers and the dense limit dominated by a condensate (see Supplementary Note for details).
We start by considering two magic-number systems composed of flexible polymers and rigid shapes, but with different dimer degeneracies: D 8R = 4 × 28 = 112 for A 8 :B 8R and D 8U = 4 × 2 = 8 for A 8 :B 8U . In the di-lute limit, we approximate the system as all dimers, in which case the combined translational entropy and internal conformational entropy of the dimers yields a freeenergy density where c is the monomer concentration of each species and L is the polymer length, and the last term captures the dependence on dimer degeneracy. In the dense limit, we treat the two species as independently adopting all possible configurations on the lattice, but then apply a correction that reflects the requirement that the occupied sites perfectly overlap, which in a mean-field approximation yields where C 8 = 2172 is the number of conformations of a flexible polymer of length L = 8 with a defined head site, while C 8R = 4 is the corresponding number of conformations for the rigid shape.
In the thermodynamic limit, the lower of the two free energies dominates. The non-convexity of the resulting free-energy density as a function of concentration implies a region of coexistence of a dilute and a dense phase, and the phase boundaries can be determined from the convex hull of the free energy ( Supplementary Fig. 3).
Varying the dimer degeneracy changes the free energy per dimer, captured by the last term in Eq. (1) and thus influences the phase boundary. The simple theory predicts a concentration of 0.61 for the phase transition of A 8 :B 8U with its dimer degeneracy of 8, while predicting that A 8 :B 8R with dimer degeneracy 112 never forms condensates. While the theory is too simplified to be quantitative, it qualitatively correctly captures the pronounced shift to higher concentration of the A 8 :B 8R phase boundary because of its higher dimer degeneracy.
The above model can be readily adapted to systems with two flexible polymers each of length L, and provides insight into the dependence of the magic-number effect on valency (see Supplementary Note). We continue to assume equal concentration of monomers and the strong interaction limit, in which case the free energy in the dilute limit is given by Eq. (1) with dimer degeneracy D L , and in the dense limit by Eq. (2), with both conformation numbers replaced by C 8 .
For the non-magic number system A n−1 :B n in the dense limit, a version of Eq. (2) that takes into account the different polymer lengths still holds (Eq. S25). In the dilute limit, the smallest fully-bonded oligomer includes n(n−1) monomers of each type. We roughly approximate the oligomer conformational degeneracy by assuming the oligomer adopts an n × (n − 1)-rectangular shape to obtain a dilute-limit free-energy density (Eq. S28). Using these free energies to compare the magic-number systems A n :B n and the non-magic-number systems A n−1 :B n , the predicted phase boundary decreases from 0.68 to 0.15 for n = 4, from 0.72 to 0.02 for n = 6, and from 0.76 to 0.01 for n = 8. Thus, this simple theory captures the trend that the phase boundary decreases more, i.e., the magicnumber effect is stronger, for higher-valency systems.

Discussion
Motivated by the key roles played by membrane-free organelles in a variety of cellular functions, we studied phase-separating systems composed of two species of multivalent polymers that form one-to-one specific bonds. In particular, motivated by the Rubisco-EPYC1 system of the algal pyrenoid, we focused on the role of rigidity of one of the multivalent components. The results reported here are based on a simple lattice polymer model supported by equally simple analytical theory; nevertheless, the model and theory capture essential features of real systems: (i ) There is a phase transition from a uniform solution to liquid droplets. (ii ) The low concentration phase is dominated by small oligomers and the high concentration phase is a gel-like condensate. (iii ) One-to-one specific interactions can lead to a magic-number effect that implies striking exceptions to the general rule that higher valency favors condensation. This last feature is a key difference with respect to standard phase-separating systems, e.g. those described by Flory-Huggins theory, in which interactions are non-specific and non-saturable, and for which there is no magic-number effect. We found that rigidity enhances the magic-number effect by increasing the relative conformational entropy of the small oligomers. Our lattice simulations are intended to provide conceptual insight, not to capture the details of polymer shapes, sizes, range of interactions, or entanglement. Nonetheless, we expect the magic-number effect and our conclusions regarding the role of rigidity to be robust with respect to these considerations, and the effect has been verified in 3D off-lattice simulations 10 .
There remain open theoretical questions. What is the nature of the condensed phase without attractive nonspecific interactions or when such interactions are repulsive? How is the phase behavior affected by more general magic-number relations, e.g. one polymer species with valency n and two partnering species with valencies p and q, with n = p + q, in particular when one or more of these species is rigid? Future work will address these questions as well as the dynamical properties of multivalent, multicomponent systems.
Consideration of phase-separating systems that rely on specific interactions has strong biological motivation. Firstly, multivalent systems with specific interactions allow for "orthogonal" phase separated droplets to form: the specific interactions holding together one class of droplets will typically not interfere with those holding together another class. Given the large number of distinct condensates now recognized within cells 21 , droplet orthogonality is a key consideration. Secondly, magicnumber effects allow for novel mechanisms of regulation. For example, the algal pyrenoid, which is a twocomponent phase-separated liquid organelle, is observed to rapidly dissolve prior to cell division. Chemical modification of the effective valency of one component to achieve a magic-number condition has been suggested as a possible driver of this rapid dissolution 10 .
It is important to ask if the magic-number effect is robust enough to be biologically relevant. In particular, protein stoichiometries are known to fluctuate within cells. In this regard, our finding that the magic-number effect persists away from equal monomer concentration for systems composed of a flexible polymer and a rigid shape provides an important new biological insight. Notably, this is exactly the case for the pyrenoid where a flexible polymer (EPYC1) interacts with a rigid enzyme complex (Rubisco).
The magic-number effect makes definite experimental predictions. First, comparable magic-number and nonmagic number systems will have very different phase boundaries. Second, suppression of phase separation by the magic-number effect is enhanced by higher valency and by higher dimer or small oligomer conformational entropy. Third, deviations from equal monomer concentration strongly shift the phase boundary (but with an exception if one species is rigid and the flexible species is in excess). We hope that the results and predictions presented here will stimulate exploration of magic-number effects in multivalent, multicomponent systems both in vitro and in vivo. the manuscript.

Model
Simulations were performed using a square grid system of 50 × 50 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 2D lattice; no more than one A-monomer and no more than one B-monomer can occupy a site.
Some polymers are considered to be "flexible" in which case any configuration of connected nearest-neighbor sites is allowed. We also consider cases where the Bpolymer (B 8R ) is a rigid 4×2 rectangle, and a "U-shaped" variant (B 8U ) described below.
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 22 . 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 two-species joint moves. Singleflexible-polymer moves are standard lattice-polymer local moves: the end-point move, the corner move, and the reptation move 10 . Single-rigid moves consist of onestep translations in the four cardinal directions and a 90-degree rotation around the center of the rigid shape. In the regime of strong specific bonds, the two species are typically held together by multiple specific bonds, which leads to dynamical freezing. This is more severe if one of the species is rigid, since the moves affect more binding sites. To enable the system to better explore configuration space, in the case that one species is rigid, we include two-species joint moves such that connected clusters of polymers move together, without breaking any specific bonds. The joint moves consist of translating a connected cluster of the two species of polymers together or rotating the whole cluster by 90-degrees around any point. To obtain thermalized ensembles, we follow a two-step sim-ulated annealing procedure: we keep k B T constant and gradually increase bond strength. We first increase 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 is increased from 0 to 11 k B T in 0.04 k B T increments, while the non-specific bond energy is kept at 0.1 k B T . Each step of annealing is simulated with 150,000-450,000 Monte-Carlo steps to ensure complete thermalization and results are averaged over 20-100 of the resulting thermalized snapshots.
Construction of the "U-shaped" A 8 :B 8U system To test how the phase diagram is affected by dimer degeneracy, we designed a variant of the A 8 :B 8R system that minimizes the number of possible dimer conformations. Specifically, we defined a "U-shaped" variant of the 4 × 2 rectangle (B 8U ), such that an A-polymer crossing the central long axis of the rectangle (except at one end) only contributes one binding energy from the two adjacent binding sites. This energetically favors A-polymer partners that follow the U-shape ( Supplementary Fig. 1).

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 valency
Valency is a key parameter in both natural and engineered two-component systems that form specific bonds [11][12][13][14] . It is therefore natural to ask how valency 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 valencies ( 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 valency 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 (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 valency is lower for the nonmagic-number system, which competes with the magic-number effect. Indeed, for the lowest-valency case, A 3 :B 4 has a similar phase boundary to A 4 :B 4 (Fig. 2a). By contrast, for the higher-valency systems, the phase boundaries in the regime of strong interactions occur at substantially higher concentrations for the magic-number systems (Fig. 2b,c). This confirms that higher valency increases the magnitude of the magic-number effect. Note that at lower interaction energies, the presence of unbonded monomers means magic-and non-magic-number 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 450,000-9,000,000 Monte-Carlo steps. 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.

Flory-Huggins theory: a brief review
Flory-Huggins theory 19 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: (S15) 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). Specifically, there are 1 M N olig (L1L2−1) M ! N olig !(M −N olig L1L2)! center-of-mass configurations, where N olig = cM/(L 1 L 2 ) is the number of oligomers. Each oligomer has a degeneracy D olig = 2D 1 D 2 , where the factor of 2 comes from the two possible orientations of the oligomer rectangle, and D 1 and D 2 are the number of possible configurations of the two species within the rectangle. Then the total number of configurations is: and the free energy is