The solubility product extends the buffering concept to heterotypic biomolecular condensates

Biomolecular condensates are formed by liquid-liquid phase separation (LLPS) of multivalent molecules. LLPS from a single ("homotypic") constituent is governed by buffering: above a threshold, free monomer concentration is clamped, with all added molecules entering the condensed phase. However, both experiment and theory demonstrate that buffering fails for the concentration dependence of multicomponent ("heterotypic") LLPS. Using network-free stochastic modeling, we demonstrate that LLPS can be described by the solubility product constant (Ksp): the product of free monomer concentrations, accounting for the ideal stoichiometries governed by the valencies, displays a threshold above which additional monomers are funneled into large clusters; this reduces to simple buffering for homotypic systems. The Ksp regulates the composition of the dilute phase for a wide range of valencies and stoichiometries. The role of Ksp is further supported by coarse-grained spatial particle simulations. Thus, the solubility product offers a general formulation for the concentration dependence of LLPS.


Introduction
Biomolecular condensates comprise a novel class of intracellular structures formed by a biophysical phenomenon called liquid-liquid phase separation (LLPS) (Banani et al., 2017;Hyman et al., 2014;. These structures serve as membraneless compartments where complex biochemistry can be organized and facilitated ; for example, T cell receptor-mediated actin nucleation efficacy spikes up multifold when all the associated signaling molecules concentrate into a condensate (Su et al., 2016) near the plasma membrane. These structures are also implicated in many age-related or neurological diseases Alberti et al., 2019;Mathieu et al., 2020).
Numerous theoretical and experimental studies have illuminated many of the biophysical requirements for condensate formation (Li et al., 2012;Wang et al., 2018). In particular, it is firmly established that clustering of weakly interacting multivalent proteins or nucleic acids is a prerequisite for the phase separation. Even a sufficiently concentrated solution of a single self-interacting protein (homotypic interaction) with multiple binding sites in its sequence can partition into protein-dense and dilute phases. Such homotypic systems display a strict threshold concentration above which phase separation occurs. This phase separation serves as a buffering mechanism for the protein in the dilute phase Klosin et al., 2020), which remains at the threshold concentration; thus, as more protein is added to the system, the dense phase droplets grow in size and number, keeping the concentration clamped in the dilute phase. Although a homotypic system closely conforms to a single fixed threshold concentration, the picture gets much more complex with multicomponent (heterotypic interactions) systems, which underlie all the biomolecular condensates found in living cells and contain complex mixtures of multivalent proteins and/or nucleic acid. Detailed theoretical analysis of lattice-based simulations explained why the dilute phase concentrations of specific components need not stay fixed when phase separation is driven by heterotypic interactions (Choi et al., 2019). This was recently followed by a thorough experimental study of the thermodynamics of the liquid-liquid phase transitions in heterotypic systems, showing clearly that concentration thresholds for phase separation no longer remain fixed and vary with relative compositions of interacting binding partners (Riback et al., 2020). But a theoretical framework for quantitatively predicting such complex varying concentration thresholds is still lacking.
We have previously distinguished between strong multivalent interactions, which can produce molecular machines, and weak multivalent interactions, which can produce what we termed 'pleomorphic ensembles' (Mayer et al., 2009;Falkenberg et al., 2013). The strong binding affinities in molecular machines (e.g., ribosomes, flagella) enforce a specific parts list with a fixed stoichiometry. Pleomorphic ensembles (e.g., cytoskeletal polymers and their associated binding proteins, neuronal post-synaptic densities, etc.) are much more plastic in their composition than molecular machines. Biomolecular condensates are a subclass of pleomorphic ensembles in that their molecular components simultaneously coexist within both a distinct phase and as solutes in the surrounding solution. In trying to quantitatively understand the relationships governing condensate formation from heterotypic components, we asked if there may be lessons to learn from ionic solution chemistry, where salts in solution are in equilibrium with a solid crystalline phase composed of a lattice of counter ions. We realized that the liquid phase separation at threshold concentrations of heterotypic biological molecules might resemble the precipitation of anions and cations from solution.
Consider a salt (let's say silver chloride, AgCl) in water; dissolution takes place until the solution becomes saturated; further addition of salt results in precipitation. In the simplest case of pure AgCl, this seems to be similar to the behavior of a homotypic (single component) condensate; that is, above a threshold total concentration of dissolved AgCl, the salt will not dissolve further, maintaining a clamped concentration of Ag + and Clin the solution above, whatever amount of AgCl is added. However, the key concept is that the saturation threshold is governed by the solubility product (SP) -the product of the individual concentrations of Ag + and Cl -, [Ag + ] * [Cl -]. Precipitation starts when this product reaches the thermodynamic parameter called the solubility productconstant, Ksp. For AgCl Ksp = 1.7 Â 10 À10 M 2 at 25˚C. Importantly, if we add any Clin the form of a highly soluble salt (e.g., KCl) into the solution, some AgCl will precipitate to maintain the solubility product at Ksp. Of course, the crystal lattice in an ionic solid is very different from the set of weak multivalent interactions inside a biomolecular condensate. But we wondered whether the Ksp might at least approximately be used to understand heterotypic LLPS.
We explore how well Ksp may apply to biomolecular condensates using two stochastic modeling approaches. With a non-spatial network-free simulator (NFsim; Sneddon et al., 2011), we show how Ksp can approximately predict the phase separation threshold for a two-component system by systematically changing the concentration of individual components with a variety of valencies. We simulate the dynamics of cluster formation, demonstrating a dramatic transition from a stable distribution of small oligomers below Ksp to an unstable bimodal distribution of small oligomers and explosively growing large polymers at or above Ksp. We will also show how a more complex mixedvalent three-component system also conforms to a Ksp. Moreover, these simulation results help explain experiments (Riback et al., 2020), where individual components of heterotypic biomolecular condensates are not effectively buffered in the dilute phase. We then expand on prior work on the structural features governing multivalent clustering (Chattaraj et al., 2019) with a spatial kinetic simulator (SpringSaLaD; Michalski and Loew, 2016) to support these conclusions and also point to some limitations.

Results
Phase boundary of a two-component system with molecules of the same valency Our baseline model system consists of a pair of tetravalent molecules (A 4 and B 4 ), each having four binding sites that can interact with an affinity (Kd) of 350 mM. We start the NFsim simulation with the same concentrations of monomers of each type and measure the free monomer concentrations when the system reaches the steady state ( Figure 1A). The product of both free concentrations is called solubility product (SP). As we increase the total concentration (synchronously of both A 4 and B 4 ), the concentration of monomeric molecules initially goes up ( Figure 1A, inset) and of course SP also goes up; upon reaching a concentration threshold (total concentration~120 mM in this case), SP plateaus to a constant value (169 mM 2 ). This is the solubility product constant or Ksp for this pair of tetravalent molecules.
Next, we ask how the system might behave differently below and above the Ksp. Similar to an analysis we employed in previous work using spatial simulations (Chattaraj et al., 2019), we probed for how the number of available molecules (i.e., the system volume) affects the cluster size distribution below and above the threshold, 80 mM and 160 mM ( Figure 1B, C). If the system shows no tendency to condense into large clusters, the size distribution will be insensitive to the number of molecules at a given total concentration. Figure 1B shows this to be the case for the 80 mM total concentration (SP = 154 mM 2 < Ksp). The shape of the cluster size distribution displays an exponential decline from monomers to higher oligomers, and this shape is insensitive to increasing the number of molecules (i.e., volume) in the system (unbinned histograms are in Figure 1-figure supplement 1A). Even in the presence of 800 molecules, there are hardly any clusters greater than 40 molecules (lowest panel of Figure 1B). Approximately 80% of the total molecules are in clusters containing less than 10 molecules, no matter how many molecules are available in the system. Extrapolating to a macroscopic system, this would be equivalent to a single soluble phase consisting of mainly monomers and small oligomers. Figure 1C and Figure 1-figure supplement 1B illustrate the cluster distribution for total concentration = 160 mM, above the threshold for constant SP (SP = 169 mM 2 = Ksp). The histogram does change shape, extending the tail to larger clusters eventually to become a bimodal distribution    as we feed more molecules into the system (i.e., as we increase the volume). With 1600 molecules in the system, more than 20% of the total molecules are in clusters larger than 100 molecules. Figure 1-figure supplement 2 shows simulations with 10,000 molecules, averaged over 100 trajectories, systematically varying the total concentrations by changing the simulation volume (as opposed to changing the number of molecules within a fixed volume in Figure 1A); the Ksp is still 169 mM 2 and bimodal distributions clearly develop above Ksp. We note that the long tail in these histograms (unbinned distribution in We hypothesize that this tendency to form increasingly larger clusters with more available molecules is a hallmark of phase separation behavior, as previously established (Chattaraj et al., 2019), and that a constant solubility product (Ksp) is a quantitative indicator to mark the threshold that underlies biomolecular condensates. Here, we use the percolation boundary, the threshold for forming large clusters, as a proxy for phase separation, realizing that they may not be completely coincident (Choi et al., 2020a;Choi et al., 2020b). Below a threshold total concentration, when SP has not reached the constant level of Ksp, the tendency to form large clusters is low and the system    would exist as a single phase (e.g., Figure 1B). Above the threshold where SP converges to the Ksp, the system tends to form very large clusters, yielding two different phases, manifest as bimodal cluster size histograms ( Figure 1C, Figure 1-figure supplement 2B, and Figure 1-figure supplement 3). The dense phase containing the larger clusters grows in size and the dilute phase concentration remains constant. Importantly, this behavior is precisely that of a buffering system, which has been proposed as one of the important biophysical functions of biomolecular condensates. The generality of this hypothesis will now be further explored through additional modeling scenarios.

Simulations with monomeric A 4 and B 4 maintained at fixed concentrations
We now further demonstrate that Ksp marks the phase transition threshold using an alternative modeling approach. In the first approach, we used a fixed total concentration (FTC) of molecules and measure the free monomer concentrations as the system reaches the steady state. In this second approach, we clamp the monomer concentrations to a constant value (clamped monomer concentration ['CMC,], Figure 2A) and allow the total concentration (free plus bound) to change over time. This is achieved by creating reactions that rapidly create and destroy monomers, such that the concentration is clamped at the ratio of these rate constants -as long as these rates are much faster than the rates of the binding reactions. The SP, in this case, is simply the product of CMC_A 4 and CMC_B 4 .
For CMC_A 4 = CMC_B 4 = 12.9 mM (SP = 166.4 mM 2 , below Ksp), total molecular concentrations rise up initially and then converge to a steady state (~56 mM) ( Figure 2B). However, going to CMC =  13 mM (SP = 169 mM 2 = Ksp), the total concentrations never reach a steady state. This phenomenon is more pronounced for a higher CMC (13.1 mM). Clearly the system is undergoing a fundamental change around SP = 169 mM 2 , identical to the threshold determined for FTC ( Figure 1). Both these modeling paradigms indicate that there is a solubility product constant (Ksp = 169 mM 2 ) beyond which the system has a much higher propensity to form larger molecular clusters, the prerequisite for phase-separated droplet formation. Another striking feature is the variability of the total concentrations at CMC = 13 mM as illustrated by the error bars around the mean counts (second panel, Figure 2B). When we look at the individual trajectories ( Figure 2C and Figure 2-figure supplements 1-3), the system shows large fluctuations and variable lag times before irreversible growth near the phase boundary ( Figure 2C, CMC = 13 mM); the total concentration explodes in some runs or fluctuates around a metastable state within the given time frame. The behavior is less stochastic away from the Ksp ( Figure 2B, C, CMC = 12.9 mM and 13.1 mM) as quantified by the relatively narrower error bars. The behavior at Ksp represents stochastic nucleation of clusters containing enough crosslinking that disassembly becomes unlikely; such larger clusters are sufficiently stable only at or above Ksp.
When we compare the outcomes from FTC and CMC methods, we see that the results are consistent with each other ( Figure 1A, Figure 2B, and Figure 2-figure supplement 4). Below the phase boundary, if we clamp the monomer concentrations to the value of free concentrations obtained from the FTC method, we recover the same steady-state total molecular concentrations (top six panels in Figure 2-figure supplement 4). However, at even slightly above the Ksp, the CMC total concentration increases with time, rather than leveling off to a higher steady value (bottom panels in

Phase transition depends on Ksp even when individual monomer concentrations are unequal
From ionic solution chemistry, we know that irrespective of the individual ionic concentrations, if the product of ion concentrations exceeds the Ksp of the salt (i.e., supersaturation), we always get precipitation to restore the solution concentrations to Ksp, even if one ion is present at a different concentration than the other ('common ion effect'). We wanted to test whether that simple chemical principle works for our relatively complex multivalent molecular clustering system. In Figure 3, we analyze two cases when the product of reactant's concentrations (SPs) is above (SP = 172 mM 2 ) and below (160 mM 2 ) the Ksp (169 mM 2 , derived from Figure 1). For each of those SPs, we vary the CMCs of A 4 and B 4 in such a way that the products of the CMCs are always equal to the assigned SP. Satisfyingly, we find that for SP < Ksp ( Figure 3B), the systems converge to stable steady states (no phase transition); but for SP > Ksp ( Figure 3A), irrespective of individual clamped concentrations, systems always show unbounded growth (phase transition). We also titrated unequal FTCs, maintained at a ratio of 3:2 (Figure 3-figure supplement 1). Interestingly, in this computational experiment the free monomer concentration of the lower abundant component (B 4 ) gets exhausted disproportionately as the threshold is approached, so that SP cannot quite reach Ksp (169 mM 2 ) and actually begins to diminish somewhat at still higher FTC. However, cluster size distribution (Figure 3-figure supplement 1B) still becomes increasingly bimodal with higher concentrations suggesting a phase separating behavior. Thus, even when monomer supply becomes depleted, the Ksp still serves as an upper limit for free monomer concentrations.

Higher valency promotes phase transition by reducing the Ksp
Increasing valency is known to increase the propensity for phase separation Mathieu et al., 2020). Therefore, we ask how the valency of the interacting heterotypic monomers affects the Ksp. We altered the molecular valencies (number of binding sites per molecule) from 3 to 5 and compute the SP profiles as a function of total concentrations (Figure 4). The total concentration needed to reach the Ksp goes down with higher valency, consistent with experiment. Going from 3v,3v pair to 4v,4v pair, Ksp changes over fivefold (852 mM 2 to 169 mM 2 ), whereas approximately threefold change (169 mM 2 to 55 mM 2 ) can be observed for transitioning into 5v,5v pair from 4v,4v.

Ksp for a mixed-valent system
We next explore what happens when we mix molecules with different valencies. Consider a pentaand trivalent (A 5 -B 3 ) molecular pair ( Figure 5). To optimize the clustering such that all sites could potentially be bound would require a stoichiometry of 3A 5 :5B 3 . Maintaining this concentration ratio, as we titrate up the total concentrations, we see an interesting pattern ( Figure 5A, inset): the free monomer concentration of B 3 , which is present in excess, goes up steadily; but the free A 5 goes up first and then starts to go down. When we take the product of free monomer concentrations (Figure 5-figure supplement 1), we do not see an SP profile that plateaus to a constant Ksp (as in Figure 1A). However, when we correct the SP expression by taking the ideal stoichiometry into account, SP = (free A 5 ) 3 (free B 3 ) 5 , we get an SP profile that does plateau to a fixed Ksp beyond the total concentration threshold of~128 mM ( Figure 5A). The Ksp expression for this mixed-valent binary system is analogous to a mixed-valent salt like Al 2 (SO 4 ) 3 . Indeed, examining the cluster size distribution confirms that this mixed-valent system has a concentration threshold at the same total concentration 128 mM (48 mM A 5 + 80 mM B 3 ) where this stoichiometry-adjusted SP becomes constant; beyond that point the cluster size distribution becomes bimodal and more and more molecules populate the larger clusters ( Figure 5B). Importantly, the free monomer concentrations do not display buffering above this threshold ( Figure 5A, inset), with the concentration of the pentavalent monomer actually decreasing. Thus, the Ksp analogy between ionic solution chemistry and biomolecular condensates seems to hold for even these more complex stoichiometries.

A ternary heterotypic system lacks dilute phase buffering while still being governed by Ksp
In some elegant recent experiments, Riback et al., 2020 titrated up the concentration of one component of several cellular multicomponent biomolecular condensates and showed that the expected buffering behavior did not pertain to heterotypic systems. We decided to see if our simple non- spatial simulations, which only consider binding and valency, could recapitulate these experimental observations.
To do this, we consider a mixed-valent three-component system. Component A 3 has three sites that can bind to a single domain of component B 1,3 ; six sites in component C 6 interact with three different sites on component B 1,3 . In our simulations, all these bindings are assumed to have the same weak affinities (Kd = 350 mM). We started by establishing conditions where B 1,3 and C 6 alone could form a bimodal cluster distribution ( Figure 6-figure supplement 1), corresponding to a phase separation. We found that this binary system has a Ksp of~2700 mM 3 (the units correspond to the ideal stoichiometry of 2 B 1,3 :1 C 6 ). We then chose total concentrations of 120 mM B 1,3 and 60 mM C 6 (well above the phase transition in Figure 6-figure supplement 1B) and performed a series of simulations with increasing levels of A 3 (Figure 6). Figure 6A-C shows three ways to analyze these data, which we chose to match the way experimental data for titration of NPM1, a key component of the nucleolus, was analyzed in Riback et al., 2020 (shown in the corresponding insets in the panels of Figure 6A-C). We do not know the valencies and affinities for the components that make up the nucleolus, an archetypal biomolecular condensate, so we made no attempt to match the data quantitatively. However, we are gratified with the obvious qualitative match to the experimental patterns, especially the ability of our simple binding model to show how A 3 does not display simple buffering in this scenario. Specifically, buffering, as observed in homotypic biomolecular condensates, would result in a plateau level of monomeric A 3 (or NPM1) as a function of total A 3 (the illustrative  'expected' behavior is shown for NPM1 in red in Figure 6A, inset). Instead, our simulations and the corresponding data on the levels of NPM1 in the nucleoplasm (dilute phase) show an increasing concentration of the titrant. Because the dilute phase would contain small oligomers, not just monomer, we confirmed that the patterns in Figure 6A, B for the monomeric A 3 are also present for small oligomers of A 3 (Figure 6-figure supplement 2).
Importantly, Figure 6D demonstrates how this complex ternary system obeys the Ksp just as well as the previously analyzed binary systems. The SP for this system is calculated based on the ideal valency matching stoichiometry: SP = [A 3 ] 2 [B 1,3 ] 6 [C 6 ] 3 . The blue triangles show this analysis for the same simulations that generated Figure 6A-C, where the total concentrations of B 1,3 and C 6 are kept fixed sufficiently high to be phase separated without any A 3 (Figure 6-figure supplement 1B) at 120 mM and 60 mM, respectively, while the total concentration of A 3 is titrated from 1 mM up to 200 mM. The red stars show a computational experiment where the total concentrations of the three Figure 6. Titration of one molecular component in a heterotypic condensate yields correlations with the experimental results [Riback et al., 2020]. We begin with 120 mM B 1,3 and 60 mM C 6 , which displays a bimodal cluster distribution (condensate formation; Supp.    components are varied concertedly, maintaining the ideal stoichiometric ratios of 2 A 3 :6 B 1,3 :3 C 6 . Both of these titrations plateau to the same value of Ksp of~10 12 mM 11 , even though they reach this Ksp at different values of the total concentration. The red star simulations display the characteristic bimodal cluster size distributions at a total concentration of [A 3 ] + [B 1,3 ] + [C 6 ] = (30 + 90 + 45) mM = 165 mM, that is, just when SP plateaus to the Ksp (Figure 6-figure supplement 3). Thus, even for this more complex mixed-valent ternary system, a threshold for monomer concentrations is set by a Ksp, and it determines the point for formation of the large cluster phase. Figure 6-figure supplement 3A also shows how the individual monomer concentrations (insets) continue to change dramatically even after the Ksp is reached. Taken together, these results demonstrate the usefulness of the Ksp concept in explaining the concentrations of individual components of complex heterotypic multivalent binding systems and how they lead to LLPS.

Spatial simulations
Because of the efficiency of the NFSim non-spatial stochastic simulator, we were able to rapidly explore many scenarios using large numbers of molecules and demonstrate that the solubility product constant (Ksp) may generally serve as a quantitative indicator for phase transitions of multivalent heterotypic binders. These simulations also allowed us to focus on only the effects of binding valency and stoichiometry. We now apply a spatial simulation framework, SpringSaLaD (Michalski and Loew, 2016), where the roles of spatial features, such as steric hindrance, molecular flexibility and proximity, may also impact Ksp. A biomolecule is modeled as a collection of spherical sites connected by spring-like linkers. The spheres may be designated as binding sites within a rule-based modeling interface, assigning them macroscopic on and off rates that the software translates to reaction probabilities as the spheres penetrate a computed reaction radius dependent on the onrate constant. Thus, the software is amenable to computational experiments where the geometric and reaction parameters of the system can be systematically varied.
Our spatial system consists of a pair of matched tetravalent molecules (A 4a and B 4b ) where each molecule contains four binding sites, with interspersed pairs of linker sites ( Figure 7A). These linker sites impart flexibility to the molecules mimicking the intrinsically disordered linker sequences that are found in many phase-separating multivalent proteins (Posey et al., 2018). Each of the magenta sites can bind to each of the green sites with an affinity of 350 mM. We begin with 100 A 4a and 100 B 4b molecules, randomly distributed in a three-dimensional rectangular volume. As the system relaxes to the steady state, we quantify the monomer concentrations and plot their product (SP) as a function of the initial total concentration -the FTC approach we described above for NFsim. We generate an SP profile by systematically varying the FTC by changing the volume of the compartment, keeping the total molecular numbers same ( Figure 7B). The SP converges to a constant value (Ksp = 318 mM 2 ) at a threshold FTC (~138 mM in this case). When we look at the detailed cluster size distributions at steady state (Figure 7-figure supplement 1A), the dimer emerges as a preferred configuration both below and above the Ksp. This dimer preference arises from the matching valency and spatial arrangement of the binding sites in the two partner molecules; such an effect, which we term a 'dimer trap,' cannot be realized in non-spatial methods, such as NFsim. Above Ksp, however, most of the individual trajectories, much like NFsim, produce discontinuous distributions with small oligomers along with one or two very larger clusters (Figure 7-figure supplement 2). This obvious bifurcation is obscured when averaging over multiple individual histograms (Figure 7figure supplement 1A) because the consistent population of small oligomers is reinforced while individual large clusters will become smeared out.
Having established the Ksp with the FTC spatial simulations, we turned to the CMC approach, again using high values for k create and k decay (as in Figure 2A). One can think of the CMC approach as being equivalent to having a large external reservoir of monomers that can rapidly diffuse into a volume where clustering is enabled. As we titrate up the CMCs (Figure 7C), the system undergoes a fundamental change at 17.85 mM, corresponding to SP = 318.6 mM 2 , which is approximately at the Ksp we derived from the FTC calculations. Below that boundary (CMC = 16.7 mM), total concentrations converge to a steady state ( Figure 7C, top panel); at (CMC = 17.85 mM) or above the boundary (CMC = 18 mM), the total concentration of molecules keeps on going up with time ( Figure 7C, second and third panels). When we look at the individual total concentration trajectories for CMC = 17.85 mM (Figure 7-figure supplement 3), much like our NFsim results, some trajectories fluctuate around a lower steady state (dilute phase) for the given time frame while some trajectories shoot up     (dense phase) after a variable lag. This stochastic behavior is also present for CMC = 18 mM, although the rate of growth is much faster in this case. As quantified by the error bars ( Figure 7C, second and third panels), stochastic fluctuations are somewhat greater in spatial simulations than the non-spatial scenario. This stochasticity can be attributed to the variable lag before the nucleation of a sufficiently large cluster for irreversible and accelerating growth; higher variability in spatial simulations is expected due to a larger number of contributing factors like steric crowding and optimal geometric conformations of binding sites.
Certain features of these results are better appreciated via Figure 7-videos 1-3, each corresponding to a typical single trajectory at the three CMCs. Interactive 3D visualizations of these three simulations are available on the 'Simularium' website hosted by the Allen Institute for Cell Science; readers may access them here: Below Ksp; At Ksp; Above Ksp. The videos each show the actual dynamics of cluster formation and diffusion within the 3D volume, synchronized with the time course of total molecular concentration and a dynamic histogram of cluster size distribution. Figure 7video 1, where monomer concentration is clamped at 16.7 mM (SP = 279 mM 2 , below Ksp), shows that the total molecular concentration fluctuates around~80 mM, matching the third data point in Figure 7B. While this single trajectory is noisy, it corresponds well to the average of 50 trajectories in the top panel of Figure 7C. Importantly, the dynamic histogram in Figure 7-video 1 shows that the system rarely samples a cluster size greater than 15 molecules, which we associate with a single dilute phase. Figure 7-video 2 displays a typical trajectory with CMC at the Ksp (monomer concentration 17.85 mM, SP = Ksp = 318.6 mM 2 ). Instead of the steady state attained below Ksp (Figure 7video 1), Figure 7-video 2 displays a noisy but accelerating increase in total concentration to a maximum of 500 mM at 400 ms and a corresponding filling of the simulation volume; the dynamic histogram (note the logarithmic scale of the x-axis compared to Figure 7-video 1) shows primarily small clusters until about 200 ms, followed by a steady siphoning of newly appearing monomers into a single large cluster, which we associate with phase separation. Figure 7-video 3, with CMC set just above Ksp at 18 mM, illustrates how one trajectory reaches a metastable steady state that lasts until~180 ms, but ultimately explodes to almost 800 mM total concentration, virtually filling the available volume (as also noted for the corresponding averaged trajectories in the lowest panel of Figure 7C). This corresponds to a nucleation step, which lasts until the formation of a sufficiently large cluster to capture most newly appearing monomers. Motions and spatial locations of individual clusters can be visualized through Figure 7-video 4, which is based on the same simulations used to produce Figure 7-video 1 and Figure 7-video 3 (i.e., below and above Ksp, respectively). Figure 7-video 4 displays the individual clusters in a given time frame by computing their centroids and a radius of gyration around that center. For visual clarity, the cluster sizes are scaled down proportionately (by a factor of 4); for example, the largest cluster in the last time frame of the above Ksp (right panel) has a radius of gyration of~48 nm. The dimension of the simulation volume is 100 * 100 * 120 nm 3 . Through Figure 7-video 4, we can better appreciate the dramatically different spatial distributions of clusters below and above Ksp. The left panel, below Ksp, corresponds to a collection of small clusters homogeneously distributed across the simulation volume (single dilute phase), while the right panel illustrates the evolution of a large cluster (dense phase) coexisting with a pool of randomly distributed small clusters (dilute phase).
Like our non-spatial simulations, both FTC and CMC approaches yield self-consistent results for the spatial system: as the monomer SP remains below the Ksp threshold (318 mM 2 ), the system exhibits only one phase; but above that threshold boundary, the molecules get partitioned into two different phases -a dilute phase with monomers and small oligomers and a highly clustered phase ( Figure 7D and Figure 7-figure supplement 1B). The validity of the Ksp is preserved even when the CMCs are unequal. We illustrate this by choosing two SPs, 280 mM 2 and 330 mM 2 , respectively below and above the Ksp, and varying the individual CMCs ( Figure 7E, F). Both the cases with SP = 280 mM 2 converge to a steady state, while SP = 330 mM 2 combinations explode in both cases. It should be noted that as the total concentration explodes in these simulations, the volume can become filled with newly created molecules; this puts a brake on the total concentration in long duration simulations.

Discussion
It has been widely understood that for single component, self-interacting (homotypic) multivalent systems undergoing liquid-liquid phase separation (LLPS), the concentration of monomer in the dilute phase remains fixed above the phase transition, no matter how much of the monomer is added to the system; this feature of biomolecular condensates underlies buffering and noise reduction Klosin et al., 2020). Of course, the maintenance of a fixed concentration in a saturated solution of a solute in equilibrium with its solid phase is an elementary thermodynamic rule. Recognizing this, we set out to see how far the analogy to solution chemistry could take us in considering heterotypic interactions between different multivalent binders. Somewhat more complex than the saturated solution of a single solute is the precipitation of solid salts from saturated solutions of their ions. Chemists well know that ionic solutions of weakly soluble salts are governed by the solubility product constant, Ksp = [C m+ ] n [A n-] m , where m and n are the valencies, respectively, of the cation C m+ and anion A n-; importantly, n and m are the stoichiometries, respectively, of the cation and the anion in the solid phase to balance the positive to the negative charges. The solubility product constant derives from a fundamental thermodynamic principleequality of chemical potential between coexisting phases (ions in solution and solid). We wondered whether similar expressions could define the thresholds for LLPS in multicomponent (heterotypic) multivalent interaction systems.
One limitation in the analogy is that the composition of the ionic solid phase, and therefore the activity, is invariant, making the system-free energy dependent on only the activities of the ions in solution. However, the ideal stoichiometry, which absolutely constrains the composition of an ionic solid, is not so strictly enforced in the condensed phase of a multivalent condensate because of the weak binding affinities that underlie these systems. Therefore, to explore how well the Ksp might apply to LLPS, we used a non-spatial stochastic network fee simulator, NFsim (Sneddon et al., 2011); it isolates only the effect of valency and binding on the concentration dependence of clustering. We used a single weak binding affinity (Kd = 350 mM) for all our simulations. The efficiency of this computational method also made it possible to screen many scenarios with a sufficient number of molecules and trajectories to assure statistically that we were not missing any interesting effects. We used two approaches to assess the threshold behavior. First, we titrated up the total concentration of pairs of multivalent binders. We found that there was a threshold above which the concentration of free monomers obeyed a Ksp expression -that is, [A] n [B] m , where m and n are the ideal stoichiometries for an oligomer with all binding sites occupied. Below the Ksp threshold, the histogram of cluster size distributions tails off exponentially from monomers to small oligomers and its shape is independent of the number of molecules available (e.g., Figure 1B and Figure 1-figure supplement 1A). However, for total concentrations above the Ksp threshold, the histogram of cluster size distributions becomes bimodal, with an increasing population of huge clusters as more molecules are added to the system (e.g., Figure 1C and Figure 1-figure supplement 1B); we consider this bimodal cluster size distribution to be a hallmark of phase separation. The individual trajectories show a separation between small oligomers and a single large cluster (Figure 1-figure supplement  3). This bifurcation between one large cluster and small oligomers has been used to define a percolation boundary, generally considered a proxy for phase separation (Choi et al., 2019;Choi et al., 2020a). To further relate the behavior of these stochastic systems to the macroscopic phase transition, we used a second modeling approach where we clamped the monomer concentrations (CMC) while allowing the clusters to grow. When monomer concentrations were clamped below the Ksp, the system reached a steady state identical to that of the corresponding closed fixed total concentration simulations (Figure 2-figure supplement 4), with identical cluster distributions containing a single decaying histogram of cluster sizes. However, when the CMC was set to the Ksp or slightly above it, the plot of total concentration vs. time exploded following an initial lag (e.g., Figure 2), indicating that as new monomers enter the system they were funneled into large clusters without attaining a steady state. Together, we feel that the behavior of these two different modeling approaches shows that the Ksp defines a threshold of monomer concentrations above which a phase separation occurs.
The generality of the Ksp as an indicator of threshold was then further tested using different scenarios. We tested a situation where the total concentrations of each molecule in the binary tetravalent heterotypic system were unequal; while these simulations showed that the Ksp determined in the analysis of the equal concentration system was obeyed for a range of unequal concentrations (Figure 3), there will be deviations if the range is extended too far (Figure 3-figure supplement  1). This deviation is not surprising considering that the binding affinities are weak, leaving an excess of binding sites empty within a cluster, when the free concentration of one binding partner is too highly depleted. We tested cases where the valencies were not identical for a binary ( Figure 5) and even a ternary system of mixed valency interactors ( Figure 6D). In these systems, a more complex Ksp formulation was required to account for the appropriate ideal stoichiometry of the cluster. In both of these cases, the Ksp was successful in defining a threshold for phase separation.
The ternary system of Figure 6 also allowed us to explore how well the simple multivalent binding simulations could reproduce some recent experiments demonstrating how the components of heterotypic biomolecular condensates are not effectively buffered in the dilute phase (Riback et al., 2020). For several prototypical cellular LLPS systems, this study elegantly showed that titration of a single component produced a free concentration of that component that increased even more rapidly than its total concentration ( Figure 6A inset is an example). The simulation results in Figure 6 show that our simple multivalent binding system is able to reproduce these experimental titration patterns; importantly, it also further demonstrates, despite this failure of simple buffering, that the concentrations of individual components in the dilute phase are still constrained by Ksp.
Recent computational (Choi et al., 2019) and theoretical (Deviri and Safran, 2021) studies demonstrated that buffering of dilute phase concentrations in multicomponent systems has a complex relationship with the interplay of homotypic and heterotypic interactions. For a two-component heterotypic system, plotting the total concentrations of one component against the other produces a phase diagram with an elliptical region corresponding to the coexistence of the dilute and condensed phases; the system has a single phase anywhere outside that ellipse. In fact, the dilute and condensed phase concentrations remain constant (i.e., buffered) along tie lines that are approximately parallel to the major axis of the ellipse; buffering fails perpendicular to tie lines (Deviri and Safran, 2021). The derivation of approximate order parameters, such as a percolation boundary (Choi et al., 2019;Choi et al., 2020b), to estimate the shape of phase diagrams, could be possible with our approach, but it is beyond the scope of this work. However, the SP should be approximately constant (Ksp) within the two-phase elliptic region. That is, the SP may be used to predict concentrations in the dilute phase even for short traversals perpendicular to tie lines in the phase diagram.
We turned to coarse-grained spatial simulations of clustering to determine if the Ksp might still be generally applicable when the shape and flexibility of multivalent molecules is explicitly considered. We used SpringSaLaD (Michalski and Loew, 2016) software; it models molecules as a series of linked spheres to represent domains within macromolecules. We had previously used this software to address the structural features that control clustering of multivalent molecules and showed that above a concentration threshold, the system display unlimited growth characteristic of LLPS (Chattaraj et al., 2019). In the present study, we examined a heterotypic pair of tetravalent binders (Figure 7), similar to the non-spatial model used in Figures 1 and 2. The SpringSaLaD structures included four binding spheres, each separated by a pair of linker spheres; the distances between sites were identical within each binding partner. The results (Figure 7) show that this model displays the same behavior noted for the non-spatial model. The product of concentrations of monomer becomes independent of the total concentration above a threshold ( Figure 7B). When monomer concentrations are clamped at or above this Ksp threshold, the total concentration explodes after a lag time ( Figure 7C, F) and the histogram of cluster sizes becomes bimodal ( Figure 7E). These results are dramatically displayed in Figure 7-videos 1-3, corresponding to typical single trajectories for the clamped monomer concentrations below, at, and above Ksp, respectively; Figure 7video 4 offers a view directly comparing cluster sizes as a function of time below and above Ksp. Thus, the Ksp concept remains valid for the spatial system shown in Figure 7A.
Manipulation of computational models thus allowed us to systematically determine how well the Ksp concept, familiar from ionic solution chemistry, might apply to the very different situation of weak interactions between multivalent macromolecules. We found that for most scenarios, if not all, the Ksp does define a threshold for the unbounded growth of large clusters and a quantitative metric for the tendency of a system to phase separate. Additionally, for those cases where there are deviations from the stereotypical behavior, insights may be gleaned into the underlying molecular factors inhibiting cluster growth. But it is important to appreciate that experimental systems may introduce additional levels of complexity, such as nonspecific low-affinity binding interactions and long-range electrostatic forces. For real biomolecular condensates, the valency and stoichiometry of the components may be unknown, and they may have multiple binding interaction of differing affinities. However, it should be possible to measure Ksp experimentally using fluorescently labeled binding partners and determining their concentrations in the dilute and condensed phase via quantitative confocal microscopy (Fink et al., 1998). Titration experiments can readily be performed in vitro and have been among the earliest studies of biomolecular condensates (Li et al., 2012). More recently, methods have been developed to perform titrations in cells (Riback et al., 2020). Such experiments could serve, initially, to validate the Ksp concept. Beyond validation, enough such data may make it possible to find optimal fits to a Ksp expression, providing an indication of the effective valency of the interactions in complex multivalent biomolecular condensates.

Non-spatial simulations (NFsim)
To develop models that probe for the effects of valency and concentrations but do not account for spatial effects, we employ the NFsim (Sneddon et al., 2011) -a non-spatial rule-based stochastic simulation framework where each biomolecule represents a molecular object that may have multiple binding sites. These sites can bind with other sites depending on a set of rules defined in the model. The simulations have units of molecular counts, but these can be readily converted to equivalent concentrations, which is how we present our results.
The NFsim model file is specified in BioNetGen Language (BNGL; http://bionetgen.org/). Let us take the example of tetravalent binders -A 4 (a1, a2, a3, a4) and B 4 (b1, b2, b3, b4). We need 16 binding rules to define all the bimolecular interactions, each having an affinity (Kd) of 350 mM. Now 1000 molecules of A 4 and B 4 with an affinity of 3500 molecules would translate to 100 mM molecular concentrations with 350 mM binding affinity. There are two equivalent ways to change the molecular concentrations: (1) change the Kd, keeping the molecular counts same, which is mathematically equivalent to changing the volume of the system; (2) change the molecular counts, keeping the Kd same. We utilize both approaches for our simulations and specify which is used in our descriptions of the results. We chose binding rules to only allow inter-molecular binding; we felt this was appropriate because NFsim cannot account for spatial proximity of binding sites or steric crowding within clusters. Once the BNGL file is defined, we then generate the corresponding XML file, which serves as the NFsim input file. We run multiple stochastic simulations in parallel using the high-performance computing facility at UConn Health (https://health.uconn.edu/high-performance-computing/). A single NFsim run (500 ms, FTC approach), containing 1000 A 4 and B 4 molecules each, took~1 min with 100 trajectories run in parallel. A Python script is then used to perform statistical analysis across all the trajectories.

Spatial simulations (SpringSaLaD)
To account for realistic spatial geometry, we employ SringSaLaD (Michalski and Loew, 2016) -a particle-based spatial simulation platform where each biomolecule is modeled as a collection of hard spheres connected by stiff spring-like linkers. The simulation algorithms are fully described (Michalski and Loew, 2016), and we previously studied various spatial biophysical factors in the context of multivalent biomolecular cluster formation (Chattaraj et al., 2019) with this software. Spring-SaLaD also uses a rule-based method to define binding reactions between multivalent binders.
The SpringSaLaD model files are generated using the graphical user interface (GUI) of the software (https://vcell.org/ssalad). We define the size of binding sites, distance between the binding sites and the overall shape of the molecule inside the GUI. To build a spatial version of the reference system (A 4a and B 4b ), two linear tetravalent molecules are constructed first, each having four binding sites and six linker sites. Unlike NFsim, in SpringSaLaD, one binding rule between 'A_type' and 'B_type' sites can take care of all the possible binding interactions. Also, we can define the binding affinity in concentration units (350 mM) directly inside the GUI. We initialize our system in a 3D rectangular geometry; for example, 100 molecules in a 10 6 nm 3 cubic volume (X = Y = Z = 100 nm) would correspond to 166 mM. We change the volume of our system to alter the molecular concentrations. Once the model is specified, as before, we run multiple stochastic simulations in parallel using our high-performance computing facility. Execution time is very much sensitive to the number of total sites due to the computational overhead of tracking individual site locations. A typical run (50 ms, FTC approach), containing 100 molecules each of A 4a and B 4b (total sites = 2000), took~6 hr.

Data analysis and visualization
Python scripts are used to analyze and visualize the data. All the scripts are written with Spyder IDE (version 4.0.0) (https://www.spyder-ide.org/). Frequently used Python libraries are numpy 1.17.3, pandas 0.25.3, and matplotlib 3.1.2. All the packages are managed by anaconda package distributions (https://www.anaconda.com/).