1. Cell Biology
  2. Structural Biology and Molecular Biophysics
Download icon

The solubility product extends the buffering concept to heterotypic biomolecular condensates

  1. Aniruddha Chattaraj
  2. Michael L Blinov
  3. Leslie M Loew  Is a corresponding author
  1. R. D. Berlin Center for Cell Analysis and Modeling, University of Connecticut School of Medicine, United States
Research Article
  • Cited 0
  • Views 1,020
  • Annotations
Cite this article as: eLife 2021;10:e67176 doi: 10.7554/eLife.67176

Abstract

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; Shin and Brangwynne, 2017). These structures serve as membraneless compartments where complex biochemistry can be organized and facilitated (Holehouse and Pappu, 2018); 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 (Shin and Brangwynne, 2017; 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; Shin et al., 2017; 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 (Holehouse and Pappu, 2018; 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 Cl- in 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 M2 at 25°C. Importantly, if we add any Cl- in 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 mixed-valent 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 (A4 and B4), each having four binding sites that can interact with an affinity (Kd) of 350 µM. 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 A4 and B4), 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 µM in this case), SP plateaus to a constant value (169 µM2). This is the solubility product constant or Ksp for this pair of tetravalent molecules.

Figure 1 with 3 supplements see all
The solubility product constant corresponds to a threshold above which molecules distribute into large clusters.

These simulation results correspond to equal total concentrations of a heterotypic tetravalent pair of molecules with Kd for individual binding of 350 μM. (A) Product of the free monomer concentrations (solubility product) as a function of the total molecular concentrations. The black dashed line indicates the plateau, corresponding to the solubility product constant (Ksp), 169 µM2. Inset plot shows the change of free molecular concentrations of both tetravalent molecules with their respective total concentrations. Each data point is an average of steady-state values from 200 trajectories. In these simulations, we titrate up the molecular counts (200, 400, 600, ...., 2000 molecules, respectively), keeping the system’s volume fixed. (B, C) Distribution of cluster sizes with varying system sizes at two different total concentrations, 80 µM and 160 µM, respectively below and above the plateau in (A). The histograms show how the molecules are distributed across different ranges of cluster sizes.

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 µM and 160 µM (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 µM total concentration (SP = 154 µM2 < 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 µM, above the threshold for constant SP (SP = 169 µM2 = 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 µM2 and bimodal distributions clearly develop above Ksp. We note that the long tail in these histograms (unbinned distribution in Figure 1—figure supplement 1B) is an average of 100 stochastic trajectories; examination of individual histograms in Figure 1—figure supplement 3 (note the logarithmic scale on the abscissa) shows almost all of them individually containing just one huge cluster along with small oligomers. Thus, Figure 1C, Figure 1—figure supplement 2B, and Figure 1—figure supplement 3 demonstrate that if the system is above Ksp, molecules are funneled into macroclusters.

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 A4 and B4 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_A4 and CMC_B4.

Figure 2 with 4 supplements see all
An alternative approach to quantify the phase transition boundary.

(A) Illustration of the clamped monomer concentration (CMC) approach. Both the molecules (A4 in magenta and B4 in green) can enter the simulation box with a rate constant kcreate (molecules/s) and exit with a rate constant kdecay (s-1). The ratio of these two parameters clamps the monomer concentration to (kcreate/kdecay). (B) Average time course (over 100 trajectories) of total molecular concentrations as a function of different CMCs. Error bars show the standard deviations across 100 trajectories. (C) Eight sample trajectories for different CMCs.

For CMC_A4 = CMC_B4 = 12.9 µM (SP = 166.4 µM2, below Ksp), total molecular concentrations rise up initially and then converge to a steady state (~56 µM) (Figure 2B). However, going to CMC = 13 µM (SP = 169 µM2 = Ksp), the total concentrations never reach a steady state. This phenomenon is more pronounced for a higher CMC (13.1 µM). Clearly the system is undergoing a fundamental change around SP = 169 µM2, identical to the threshold determined for FTC (Figure 1). Both these modeling paradigms indicate that there is a solubility product constant (Ksp = 169 µM2) 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 µM 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 13), the system shows large fluctuations and variable lag times before irreversible growth near the phase boundary (Figure 2C, CMC = 13 µM); 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 µM and 13.1 µM) 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 Figure 2—figure supplement 4).

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 µM2) and below (160 µM2) the Ksp (169 µM2, derived from Figure 1). For each of those SPs, we vary the CMCs of A4 and B4 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 (B4) gets exhausted disproportionately as the threshold is approached, so that SP cannot quite reach Ksp (169 µM2) 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.

Figure 3 with 1 supplement see all
The Ksp defines a threshold for unlimited growth of clusters even when the individual concentrations of heterotypic multivalent binding partners are unequal.

With Ksp determined from Figures 1 and 2 at ~169 μM2, solubility product (SP) was clamped above in (A) at 172 μM2 and below in (B) at 160 μM2. The solid lines (magenta and green) and error bars represent the mean and standard deviations across 100 trajectories.

Higher valency promotes phase transition by reducing the Ksp

Increasing valency is known to increase the propensity for phase separation (Holehouse and Pappu, 2018; 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 µM2 to 169 µM2), whereas approximately threefold change (169 µM2 to 55 µM2) can be observed for transitioning into 5v,5v pair from 4v,4v.

Change of Ksp with molecular valency.

For 3,3 case (blue stars), molecular counts of each type = [500, 1000, 1500, ..., 3000]. For 4,4 case (red stars), molecular counts of each type = [100, 200, 300, ..., 1000]. For 5,5 case (cyan stars), molecular counts of each type = [50, 100, 150, ..., 500]. Kd is set to 3500 molecules in all these cases. Horizontal dashed lines indicate the Ksp of the corresponding system.

Ksp for a mixed-valent system

We next explore what happens when we mix molecules with different valencies. Consider a penta- and trivalent (A5–B3) molecular pair (Figure 5). To optimize the clustering such that all sites could potentially be bound would require a stoichiometry of 3A5:5B3. Maintaining this concentration ratio, as we titrate up the total concentrations, we see an interesting pattern (Figure 5A, inset): the free monomer concentration of B3, which is present in excess, goes up steadily; but the free A5 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 A5)3(free B3)5, we get an SP profile that does plateau to a fixed Ksp beyond the total concentration threshold of ~128 μM (Figure 5A). The Ksp expression for this mixed-valent binary system is analogous to a mixed-valent salt like Al2(SO4)3. Indeed, examining the cluster size distribution confirms that this mixed-valent system has a concentration threshold at the same total concentration 128 µM (48 µM A5 + 80 µM B3) 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.

Figure 5 with 1 supplement see all
A mixed-valent binary system obeys a stoichiometry-adjusted Ksp.

(A) Logarithm of solubility product (SP) as a function of total concentrations (A5 + B3). Inset shows the variation of free molecular concentrations w.r.t. their initial total concentrations. Molecules are added at a fixed volume to vary the total concentrations. (B) Cluster size distributions become more bimodal as we go beyond the critical concentration (128 µM in this case).

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 A3 has three sites that can bind to a single domain of component B1,3; six sites in component C6 interact with three different sites on component B1,3. In our simulations, all these bindings are assumed to have the same weak affinities (Kd = 350 µM). We started by establishing conditions where B1,3 and C6 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 µM3 (the units correspond to the ideal stoichiometry of 2 B1,3:1 C6). We then chose total concentrations of 120 µM B1,3 and 60 µM C6 (well above the phase transition in Figure 6—figure supplement 1B) and performed a series of simulations with increasing levels of A3 (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 Adoes not display simple buffering in this scenario. Specifically, buffering, as observed in homotypic biomolecular condensates, would result in a plateau level of monomeric A3 (or NPM1) as a function of total A3 (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 A3 are also present for small oligomers of A3 (Figure 6—figure supplement 2).

Figure 6 with 3 supplements see all
Titration of one molecular component in a heterotypic condensate yields correlations with the experimental results [Riback et al., 2020].

We begin with 120 µM B1,3 and 60 µM C6, which displays a bimodal cluster distribution (condensate formation; Supp. Fig. S6); we then titrate up A3 concentration from 1 µM to 200 µM. (AC) Free A3 (monomeric) concentration, bound A3 (total A3 – free A3) concentration and their ratio as a function of total A3 concentrations. Inset figures are replotted from the data reported in Riback et al., 2020. To guide the eye, we fit their experimental data to a generic function, y = a * xn where a and n are pre-exponent and exponent factors, respectively. The red lines in the inset plots demonstrate the 'expected' trend if the condensation is purely driven by homotypic interactions. (D) Blue triangles correspond to the solubility product (SP) of the three-component system when A3 is being titrated up gradually, with fixed total [B1,3] = 120 µM and [C6] = 60 µM. Red stars indicate the scenario when we simultaneously change concentrations of all three components, keeping a concentration ratio of 2:6:3.

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 = [A3]2[B1,3]6[C6]3. The blue triangles show this analysis for the same simulations that generated Figure 6A–C, where the total concentrations of B1,3 and C6 are kept fixed sufficiently high to be phase separated without any A3 (Figure 6—figure supplement 1B) at 120 µM and 60 µM, respectively, while the total concentration of A3 is titrated from 1 µM up to 200 µM. The red stars show a computational experiment where the total concentrations of the three components are varied concertedly, maintaining the ideal stoichiometric ratios of 2 A3:6 B1,3:3 C6. Both of these titrations plateau to the same value of Ksp of ~1012 µM11, 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 [A3] + [B1,3] + [C6] = (30 + 90 + 45) µM = 165 µM, 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 on-rate 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 (A4a and B4b) 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 µM. We begin with 100 A4a and 100 B4b 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 µM2) at a threshold FTC (~138 µM 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 7—figure supplement 1A) because the consistent population of small oligomers is reinforced while individual large clusters will become smeared out.

Figure 7 with 7 supplements see all
Spatial simulations demonstrate similar phase boundary behavior as network-free simulator (NFsim).

(A) SpringSaLaD representations of a pair of tetravalent binders. A4a and B4b consist of four magenta and green spherical binding sites (radius = 1.5 nm) and six orange and pink linker sites (radius = 0.75 nm). Diffusion constants for all the sites are set to 2 µm2/s. For individual binding, dissociation constant, Kd = 350 µM (Kon = 20 µM−1.s−1, Koff = 7000 s−1). Simulation time constants, dt (step size) = 10−8 s−1 , dt_spring (spring relaxation constant) = 10−9 s−1 . (B) Solubility product (SP) profile of the spatial system. We place a total of 200 molecules (100 A4a + 100 B4b) in 3D boxes with varying volumes and quantify the monomer concentrations (free A4a and free B4b) at steady states. Each data point is an average over 100 trajectories. Solubility product constant (Ksp) = 318 µM2, the horizontal dashed line. (C) Total molecular concentration profiles for three clamped monomer concentrations (CMCs). The solid lines and error bars represent the mean and standard deviation over 50 trajectories. (D) Cluster size distributions at the last time point of CMC trajectories, that is, 1000 ms for CMC = 16.7 µM, 400 ms for CMC = 17.85 µM and 18 µM. More detailed histograms without binning are shown in Figure 7—figure supplement 1. (E, F) Irrespective of individual monomer concentrations, total molecular concentrations converge to steady states as long as the solubility product (SP = 280 µM2) < Ksp (E) and diverge with time when SP (= 330 µM2) > Ksp (F). The solid lines and error bars represent the mean and standard deviation over 50 trajectories.

Having established the Ksp with the FTC spatial simulations, we turned to the CMC approach, again using high values for kcreate and kdecay (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 µM, corresponding to SP = 318.6 µM2, which is approximately at the Ksp we derived from the FTC calculations. Below that boundary (CMC = 16.7 µM), total concentrations converge to a steady state (Figure 7C, top panel); at (CMC = 17.85 µM) or above the boundary (CMC = 18 µM), 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 µM (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 µM, 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 13, 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 KspAt KspAbove 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 7—video 1, where monomer concentration is clamped at 16.7 µM (SP = 279 µM2, below Ksp), shows that the total molecular concentration fluctuates around ~80 µM, 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 µM, SP = Ksp = 318.6 µM2). Instead of the steady state attained below Ksp (Figure 7—video 1), Figure 7—video 2 displays a noisy but accelerating increase in total concentration to a maximum of 500 µM 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 µM, illustrates how one trajectory reaches a metastable steady state that lasts until ~180 ms, but ultimately explodes to almost 800 µM 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 nm3. 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 µM2), 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 µM2 and 330 µM2, respectively below and above the Ksp, and varying the individual CMCs (Figure 7E, F). Both the cases with SP = 280 µM2 converge to a steady state, while SP = 330 µM2 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 (Holehouse and Pappu, 2018; 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 = [Cm+]n[An-]m, where m and n are the valencies, respectively, of the cation Cm+ and anion An-; 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 principle – equality 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 μM) 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 13, corresponding to typical single trajectories for the clamped monomer concentrations below, at, and above Ksp, respectively; Figure 7—video 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.

Materials and methods

Non-spatial simulations (NFsim)

Request a detailed protocol

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 – A4 (a1, a2, a3, a4) and B4 (b1, b2, b3, b4). We need 16 binding rules to define all the bimolecular interactions, each having an affinity (Kd) of 350 µM. Now 1000 molecules of A4 and B4 with an affinity of 3500 molecules would translate to 100 µM molecular concentrations with 350 µM 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 A4 and B4 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)

Request a detailed protocol

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. SpringSaLaD 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 (A4a and B4b), 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 µM) directly inside the GUI. We initialize our system in a 3D rectangular geometry; for example, 100 molecules in a 106 nm3 cubic volume (X = Y = Z = 100 nm) would correspond to 166 µM. 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 A4a and B4b (total sites = 2000), took ~6 hr.

Data analysis and visualization

Request a detailed protocol

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/).

All the model files, Python scripts, and a ‘Readme’ description of all the contents are available in a public GitHub repository: https://github.com/achattaraj/Ksp_phase_separation, (copy archived at swh:1:rev:22643ca2ed21b527ccdedbe6a99c2cfc29780df8), Chattaraj, 2021 .

Data availability

All the model files, Python scripts and a "Readme" description of all the contents are available in a public GitHub repository: https://github.com/achattaraj/Ksp_phase_separation (copy archived at https://archive.softwareheritage.org/swh:1:rev:22643ca2ed21b527ccdedbe6a99c2cfc29780df8). Also source data files are given for 7 figures that are part of the manuscript.

References

Decision letter

  1. Rohit V Pappu
    Reviewing Editor; Washington University in St Louis, United States
  2. José D Faraldo-Gómez
    Senior Editor; National Heart, Lung and Blood Institute, National Institutes of Health, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Acceptance summary:

The concept of solubility products will likely be useful as a tool for experimentalists studying the prospect of buffering via condensate formation in systems driven by heterotypic interactions.

Decision letter after peer review:

Thank you for submitting your article "Solubility product governs the concentration threshold for formation of biomolecular condensates" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by José Faraldo-Gómez as the Senior Editor. All reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential Revisions (for the authors):

1. Please refocus the title, abstract and scope to highlight the proposal/finding that the solubility product "rescues" the concept of buffering even in systems where heterotypic interactions drive phase transitions. The consensus is that the solubility product might be rather useful, if it can be readily calculated, to set expectations regarding thresholding/buffering in systems with heterotypic interactions, i.e., all biomolecular condensates.

2. Please improve the overall scholarship to place the work in the appropriate context. Please note that Choi et al., previously explained why the dilute phase concentrations of specific components need not stay fixed when phase separation is driven by heterotypic interactions. Please see Figure 12 and the accompanying discussion in https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007028.

Please also note that Deviri and Safran have introduced a physical theory for how buffering can be achieved when phase separation is driven by heterotypic interactions. Please see, please cite, and please include a suitable discussion of:

https://www.biorxiv.org/content/10.1101/2021.01.05.425486v1.

3. Please eliminate all remarks about freshman chemistry or high school chemistry.

4. Please provide a mechanistic explanation for the implications of solubility products setting thresholds.

5. The consensus is that the section on the effects of linkers detracts from the focus of the manuscript. A detailed and rigorous treatment of this subject has been published in eLife by Harmon et al.,. We recommend that the revised manuscript exclude all results beyond Figure 7. This, we believe, will sharpen the focus, highlight novel aspects, and avoid overlap with work that has come before.

6. The reviewers have raised specific concerns and weaknesses regarding the solubility product. It is a useful concept that touches base with percolation transitions, but not necessarily with the rigor of concepts underlying phase transitions. Specifically, the connections between solubility products and slopes of tie lines needs to be established. We believe that this is beyond the scope of the current manuscript, but it will be very helpful to discuss this issue in a revised Discussion section. Please also note that formal order parameters are not being used in either simulation paradigm for defining phase separation.

7. Please provide a synthesis to explain how either simulation paradigm and the calculation of solubility products can be used to analyze experimental data. In doing so, please spell out the type of experimental data that will come in handy in testing the hypothesis that the solubility product sets an upper limit in systems with heterotypic interactions, thereby enabling the field to assess where the buffering capacity lies in distinct multicomponent systems.

8. And please consider presenting analyses in terms of amplitudes of fluctuations. These will likely help in paring down the number of panels that one presents for each system.

Reviewer #1 (Recommendations for the authors):There, unfortunately, are far too many concerns to be raised with regard to the technical and theoretical aspects of this work. A lot of new jargon terms are introduced without connecting to established terminology. The work on linkers is a minimalist redux of that of Harmon et al. The "dimer trap", also explained by Harmon et al., is a redux of the work of Wingreen and colleagues. The unique insight here has to do with the finding / proposal that the solubility product sets an upper limit on the joint concentrations, and hence appears to rescue the concept of buffering in systems with heterotypic interactions. Recent work from Deviri and Safran and from at least one other lab have introduced the concept of heterotypic buffering and provided rigorous thermodynamic descriptions, in the language and formalism of slopes of tie lines, focusing on chemical and mechanical equilibria to explain the full range of buffering capacities realizable as a consequence of the interplay between homotypic and heterotypic interactions. If we remove all the parts of the current manuscript that are best described as redux, then we are left with the solubility product and the rescue of buffering as the novel insight. What is lacking is a rigorous connection to the concepts of phase separation and percolation, achievable only through the calculation of full phase diagrams or at least coexisting curves and the demonstration that the solubility product is invariant along tie lines, and it should be true for all tie lines.

Reviewer #2 (Recommendations for the authors):

– The resolution of the cluster sizes is quite low for the non-spatial simulations, showing only 4 bins, whereas it is much higher for the spatial ones. Why so low?

– The movies are helpful in observing the dynamics, but what are the Spring constants?

– Does Ksp still hold in Figure 5 when the Stoichiometry of components is not 5:3?

– In Figure 6, what is the fit function to the dots?

– Do all sites retain excluded volume, even when they are bound to other sites?

– What rate constants were used? What was the integration time-step? Are the results sensitive to these choices?

– I think mentioning freshman/high school chemistry once in the paper is enough. The repeated references in the discussion to the age at which we all should have learned this start to sound like a failure of the field. A text book reference could instead emphasize that the concept is old and well established.

Reviewer #3 (Recommendations for the authors):The authors claim that the solubility product is a simple yet useful concept to study the phase boundary, but in fact, the solubility product is a specific realization of a much more general concept of the reaction quotient (Q) and the equilibrium constant (K), which is also from freshman chemistry. Hence, their simulation data just show that at low concentrations Q is smaller than K, and after Q reaches K, it remains constant. What is the benefit of using Ksp instead of more general K?

Unfortunately, it seems to me that there are almost no new findings relevant to the field of phase separation biology. Importance of the linker flexibility and the sticker spacing in phase separation have been known for a while, and the valency effect as well. However, the technical advancement (or applications) might be useful to the computational field, so I strongly recommend this paper for publication in a more technical journal for computational work. I don't think that this manuscript meets the eLife criteria.

https://doi.org/10.7554/eLife.67176.sa1

Author response

Essential Revisions (for the authors):

1. Please refocus the title, abstract and scope to highlight the proposal/finding that the solubility product "rescues" the concept of buffering even in systems where heterotypic interactions drive phase transitions. The consensus is that the solubility product might be rather useful, if it can be readily calculated, to set expectations regarding thresholding/buffering in systems with heterotypic interactions, i.e., all biomolecular condensates.

We appreciate this advice to fully focus our paper on the idea that the Ksp can serve as a more general concept than buffering. We have eliminated the additional results (Figure 8) on how the Ksp might be influenced by structural features of the interacting molecules. Of course, the Ksp reduces to simple buffering for homotypic systems, so the last comment of the editor that our idea applies to all biomolecular condensates is also helpful. We believe we have captured this refocus in our new title and abstract:

“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 multi-component (“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.”

We have extensively revised the manuscript to focus on this message, to address the other issues raised by the reviewers and to further clarify the presentation. The major changes that address the Reviewers’ comments are highlighted in red in the revised manuscript.

2. Please improve the overall scholarship to place the work in the appropriate context. Please note that Choi et al., previously explained why the dilute phase concentrations of specific components need not stay fixed when phase separation is driven by heterotypic interactions. Please see Figure 12 and the accompanying discussion in https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007028.

Please also note that Deviri and Safran have introduced a physical theory for how buffering can be achieved when phase separation is driven by heterotypic interactions. Please see, please cite, and please include a suitable discussion of:

https://www.biorxiv.org/content/10.1101/2021.01.05.425486v1.

We became aware of the very relevant work in Choi et al. as a result of Reviewer 1’s comment (thank you) and now cite it in 3 places. We were aware of the Devri and Safran preprint, which was submitted after our original preprint submission. We have updated our paper with a citation and discussion of this work (see response to item 6).

3. Please eliminate all remarks about freshman chemistry or high school chemistry.

As chemists (2 of the authors), we find it exciting that this elementary chemistry principle could be extended to help explain the complex concentration dependence of LLPS; however, we now explain this without referring to “freshman” chemistry.

4. Please provide a mechanistic explanation for the implications of solubility products setting thresholds.

We have added additional explanations within the discussion, which also further highlight some of the limitations of the Ksp idea:

“The of solubility product constant derives from a fundamental thermodynamic principle – equality of chemical potential between coexisting phases (ions in solution and solid). 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”

See also the response to item 6.

5. The consensus is that the section on the effects of linkers detracts from the focus of the manuscript. A detailed and rigorous treatment of this subject has been published in eLife by Harmon et al.,. We recommend that the revised manuscript exclude all results beyond Figure 7. This, we believe, will sharpen the focus, highlight novel aspects, and avoid overlap with work that has come before.

We have complied and believe it has sharpened the focus.

6. The reviewers have raised specific concerns and weaknesses regarding the solubility product. It is a useful concept that touches base with percolation transitions, but not necessarily with the rigor of concepts underlying phase transitions. Specifically, the connections between solubility products and slopes of tie lines needs to be established. We believe that this is beyond the scope of the current manuscript, but it will be very helpful to discuss this issue in a revised Discussion section. Please also note that formal order parameters are not being used in either simulation paradigm for defining phase separation.

We have added the following paragraph to the Discussion (referencing both Choi and Devri):

“Recent computational [13] and theoretical [23] 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 [23]. The derivation of approximate order parameters, such as a percolation boundary [13, 21], to estimate the shape of phase diagrams, could be possible with our approach, but it is beyond the scope of this work. However, the solubility product should be approximately constant (Ksp) within the two-phase elliptic region. That is, the solubility product may be used to predict concentrations in the dilute phase even for short traversals perpendicular to tie lines in the phase diagram.”

7. Please provide a synthesis to explain how either simulation paradigm and the calculation of solubility products can be used to analyze experimental data. In doing so, please spell out the type of experimental data that will come in handy in testing the hypothesis that the solubility product sets an upper limit in systems with heterotypic interactions, thereby enabling the field to assess where the buffering capacity lies in distinct multicomponent systems.

We have expanded the concluding sentences to suggest how the Ksp concept could lead to new experimental insights:

“However, it should be possible to measure Ksp experimentally using fluorescently labelled binding partners and determining their concentrations in the dilute and condensed phase via quantitative confocal microscopy [24]. Titration experiments can readily be performed in vitro and have been among the earliest studies of biomolecular condensates [9]. More recently, methods have been developed to perform titrations in cells [14]. 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.”

8. And please consider presenting analyses in terms of amplitudes of fluctuations. These will likely help in paring down the number of panels that one presents for each system.

Done. Individual trajectories displaying varying time lags for nucleation, are now available in Supporting Figures. Because of this interesting variable nucleation time, we felt it was important to include them.

Reviewer #1 (Recommendations for the authors):

There, unfortunately, are far too many concerns to be raised with regard to the technical and theoretical aspects of this work. A lot of new jargon terms are introduced without connecting to established terminology. The work on linkers is a minimalist redux of that of Harmon et al. The "dimer trap", also explained by Harmon et al., is a redux of the work of Wingreen and colleagues. The unique insight here has to do with the finding / proposal that the solubility product sets an upper limit on the joint concentrations, and hence appears to rescue the concept of buffering in systems with heterotypic interactions. Recent work from Deviri and Safran and from at least one other lab have introduced the concept of heterotypic buffering and provided rigorous thermodynamic descriptions, in the language and formalism of slopes of tie lines, focusing on chemical and mechanical equilibria to explain the full range of buffering capacities realizable as a consequence of the interplay between homotypic and heterotypic interactions. If we remove all the parts of the current manuscript that are best described as redux, then we are left with the solubility product and the rescue of buffering as the novel insight. What is lacking is a rigorous connection to the concepts of phase separation and percolation, achievable only through the calculation of full phase diagrams or at least coexisting curves and the demonstration that the solubility product is invariant along tie lines, and it should be true for all tie lines.

As described above, we have refocused the paper on how the solubility product offers a more general framework than buffering for constraining the concentrations of multivalent molecules undergoing LLPS. We also attempted to better place our work, which is more chemistry oriented, in the context of previous work on the physics of LLPS.

Reviewer #2 (Recommendations for the authors):

– The resolution of the cluster sizes is quite low for the non-spatial simulations, showing only 4 bins, whereas it is much higher for the spatial ones. Why so low?

Corrected. Importantly, we now also show histograms for single FTC simulations at steady state (Figure 1 —figure supplement 3); these show clear bifurcations between an exponentially decreasing distribution at low cluster sizes in equilibrium with a single cluster at a very high molecular count.

– The movies are helpful in observing the dynamics, but what are the Spring constants?

We use the default value throughout: dt_spring = 10-9 s-1 (see input file in the Git repository). The algorithm only uses the stiff springs to transmit the Langevin forces between linked spheres; their actual deformation is negligible. See [19] for details.

– Does Ksp still hold in Figure 5 when the Stoichiometry of components is not 5:3?

Figure 3 —figure supplement 1 shows that when we deviate from ideal stoichiometry (for A4 – B4, it is 1:1), SP reaches a maximal level close to the Ksp derived from the ideal stoichiometry. The maximal point does represent the onset of phase transition. We did not try this for a titration with a non-ideal stoichiometry for the system in Figure 5, but expect it to behave similarly.

– In Figure 6, what is the fit function to the dots?

The data is fitted to a generic function y = a*x n where a and n are pre-exponent and exponent factors respectively. This was done just to guide the eye, as in the experimental data from Riback, et al.

– Do all sites retain excluded volume, even when they are bound to other sites?

In SpringSaLaD, all the sites retain excluded volumes whether bound or unbound. In other words, at each simulation time step, the solver checks for any physical overlap, irrespective of the binding status of the sites.

– What rate constants were used? What was the integration time-step? Are the results sensitive to these choices?

All the parameters can be found in simulation input files, available in a public github directory. https://github.com/achattaraj/Ksp_phase_separation.

Kon = 20 μM-1.s-1, Koff = 7000 s-1, dt = 10-8 s-1

The time courses are of course sensitive to these choices, but the paper focuses on comparisons in steady state behavior and these are sensitive to only Koff/Kon. Of course all these values were kept the same for all simulations.

– I think mentioning freshman/high school chemistry once in the paper is enough. The repeated references in the discussion to the age at which we all should have learned this start to sound like a failure of the field. A text book reference could instead emphasize that the concept is old and well established.

We have deleted all these comments.

Reviewer #3 (Recommendations for the authors):

The authors claim that the solubility product is a simple yet useful concept to study the phase boundary, but in fact, the solubility product is a specific realization of a much more general concept of the reaction quotient (Q) and the equilibrium constant (K), which is also from freshman chemistry. Hence, their simulation data just show that at low concentrations Q is smaller than K, and after Q reaches K, it remains constant. What is the benefit of using Ksp instead of more general K?

All of the points in our FTC titration diagrams are at equilibrium. The reaction quotient (Q) refers to a system that is out of equilibrium, allowing the chemist to predict how concentrations will change as equilibrium is approached, so we do not understand the relevance of this comment. Furthermore, we don’t see how an equilibrium constant, K for these complex systems could be formulated. Because of multivalency, we are dealing with systems containing multiple cascading reversible reactions and an infinite number of possible species and also with a system that undergoes a phase transition (or percolation transition). The key point in our treatment is that the SP, which can be determined empirically, converges to an approximate constant at the concentration threshold for LLPS (or more precisely, the percolation threshold). We are sorry that Reviewer 3 missed this key point and hope that the revised manuscript will help her/him understand the novelty of our study.

Unfortunately, it seems to me that there are almost no new findings relevant to the field of phase separation biology. Importance of the linker flexibility and the sticker spacing in phase separation have been known for a while, and the valency effect as well. However, the technical advancement (or applications) might be useful to the computational field, so I strongly recommend this paper for publication in a more technical journal for computational work. I don't think that this manuscript meets the eLife criteria.

Reference

1. Banani, S.F., et al., Biomolecular condensates: organizers of cellular biochemistry. Nature Reviews Molecular Cell Biology, 2017. 18: p. 285.

2. Hyman, A.A., C.A. Weber, and F. Jülicher, Liquid-Liquid Phase Separation in Biology. Annual Review of Cell and Developmental Biology, 2014. 30(1): p. 39-58.

3. Shin, Y. and C.P. Brangwynne, Liquid phase condensation in cell physiology and disease. Science, 2017. 357(6357).

4. Holehouse, A.S. and R.V. Pappu, Functional Implications of Intracellular Phase Transitions. Biochemistry, 2018. 57(17): p. 2415-2423.

5. Su, X., et al., Phase separation of signaling molecules promotes T cell receptor signal transduction. Science (New York, N.Y.), 2016. 352(6285): p. 595-599.

6. Shin, Y. and C.P. Brangwynne, Liquid phase condensation in cell physiology and disease. Science, 2017. 357(6357): p. eaaf4382.

7. Alberti, S., A. Gladfelter, and T. Mittag, Considerations and Challenges in Studying Liquid-Liquid Phase Separation and Biomolecular Condensates. Cell, 2019. 176(3): p. 419-434.

8. Mathieu, C., R.V. Pappu, and J.P. Taylor, Beyond aggregation: Pathological phase transitions in neurodegenerative disease. Science, 2020. 370(6512): p. 56.

9. Li, P., et al., Phase transitions in the assembly of multivalent signalling proteins. Nature, 2012. 483: p. 336.

10. Shin, Y., et al., Spatiotemporal Control of Intracellular Phase Transitions Using Light-Activated optoDroplets. Cell, 2017. 168(1): p. 159-171.e14.

11. Wang, J., et al., A Molecular Grammar Governing the Driving Forces for Phase Separation of Prion-like RNA Binding Proteins. Cell, 2018. 174(3): p. 688-699.e16.

12. Klosin, A., et al., Phase separation provides a mechanism to reduce noise in cells. Science, 2020. 367(6476): p. 464.

13. Choi, J.-M., F. Dar, and R.V. Pappu, LASSI: A lattice model for simulating phase transitions of multivalent proteins. PLOS Computational Biology, 2019. 15(10): p. e1007028.

14. Riback, J.A., et al., Composition-dependent thermodynamics of intracellular phase separation. Nature, 2020. 581(7807): p. 209-214.

15. Mayer, B.J., M.L. Blinov, and L.M. Loew, Molecular machines or pleiomorphic ensembles: signaling complexes revisited. J Biol, 2009. 8(9): p. 81.1-81.8.

16. Falkenberg, Cibele V., Michael L. Blinov, and Leslie M. Loew, Pleomorphic Ensembles: Formation of Large Clusters Composed of Weakly Interacting Multivalent Molecules. Biophys J, 2013. 105(11): p. 2451-2460.

17. Sneddon, M.W., J.R. Faeder, and T. Emonet, Efficient modeling, simulation and coarse-graining of biological complexity with NFsim. Nature Methods, 2011. 8(2): p. 177-183.

18. Chattaraj, A., M. Youngstrom, and L.M. Loew, The Interplay of Structural and Cellular Biophysics Controls Clustering of Multivalent Molecules. Biophys J, 2019. 116(3): p. 560-572.

19. Michalski, P.J. and L.M. Loew, SpringSaLaD: A Spatial, Particle-Based Biochemical Simulation Platform with Excluded Volume. Biophysical journal, 2016. 110(3): p. 523-529.

20. Choi, J.-M., A.S. Holehouse, and R.V. Pappu, Physical Principles Underlying the Complex Biology of Intracellular Phase Transitions. Annual Review of Biophysics, 2020. 49(1): p. 107-133.

21. Choi, J.M., A.A. Hyman, and R.V. Pappu, Generalized models for bond percolation transitions of associative polymers. Phys Rev E, 2020. 102(4-1): p. 042403.

22. Posey, A.E., A.S. Holehouse, and R.V. Pappu, Chapter One – Phase Separation of Intrinsically Disordered Proteins, in Methods in Enzymology, E. Rhoades, Editor. 2018, Academic Press. p. 1-30.

23. Deviri, D. and S.A. Safran, Physical theory of biological noise buffering by multi-component phase separation. bioRxiv, 2021: p. 2021.01.05.425486.

24. Fink, C., F. Morgan, and L.M. Loew, Intracellular fluorescent probe concentrations by confocal microscopy. Biophys J, 1998. 75(4): p. 1648-58.

https://doi.org/10.7554/eLife.67176.sa2

Article and author information

Author details

  1. Aniruddha Chattaraj

    R. D. Berlin Center for Cell Analysis and Modeling, University of Connecticut School of Medicine, Farmington, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7105-6621
  2. Michael L Blinov

    R. D. Berlin Center for Cell Analysis and Modeling, University of Connecticut School of Medicine, Farmington, United States
    Contribution
    Resources, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9363-9705
  3. Leslie M Loew

    R. D. Berlin Center for Cell Analysis and Modeling, University of Connecticut School of Medicine, Farmington, United States
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Validation, Investigation, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    les@volt.uchc.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1851-4646

Funding

National Institute of General Medical Sciences (R24 GM137787)

  • Leslie M Loew

National Institute of General Medical Sciences (R01 GM132859)

  • Leslie M Loew

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank our colleague Boris Slepchenko for thoughtful discussions on the clamped monomer concentration approach. We thank Blair Lyons and Graham Johnson of the Allen Institute for Cell Science for developing and hosting the Simularium interactive visualizations of our SpringSaLaD simulations. We gratefully acknowledge support of this work by NIH grants R01 GM132859 and R24 GM137787 from the National Institute for General Medical Sciences.

Senior Editor

  1. José D Faraldo-Gómez, National Heart, Lung and Blood Institute, National Institutes of Health, United States

Reviewing Editor

  1. Rohit V Pappu, Washington University in St Louis, United States

Publication history

  1. Received: February 2, 2021
  2. Accepted: July 2, 2021
  3. Accepted Manuscript published: July 8, 2021 (version 1)
  4. Version of Record published: July 19, 2021 (version 2)

Copyright

© 2021, Chattaraj et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,020
    Page views
  • 119
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cell Biology
    2. Medicine
    Richard Nakamura et al.
    Research Article

    Background: Blinding reviewers to applicant identity has been proposed to reduce bias in peer review.

    Methods: This experimental test used 1200 NIH grant applications, 400 from Black investigators, 400 matched applications from White investigators, and 400 randomly selected applications from White investigators. Applications were reviewed by mail in standard and redacted formats.

    Results: Redaction reduced, but did not eliminate, reviewers' ability to correctly guess features of identity. The primary, pre-registered analysis hypothesized a differential effect of redaction according to investigator race in the matched applications. A set of secondary analyses (not pre-registered) used the randomly selected applications from White scientists and tested the same interaction. Both analyses revealed similar effects: Standard format applications from White investigators scored better than those from Black investigators. Redaction cut the size of the difference by about half (e.g. from a Cohen's d of 0.20 to 0.10 in matched applications); redaction caused applications from White scientists to score worse but had no effect on scores for Black applications.

    Conclusions: Grant-writing considerations and halo effects are discussed as competing explanations for this pattern. The findings support further evaluation of peer review models that diminish the influence of applicant identity.

    Funding: Funding was provided by the NIH.

    1. Cell Biology
    2. Stem Cells and Regenerative Medicine
    Nicholas Strash et al.
    Research Article

    Multiple mitogenic pathways capable of promoting mammalian cardiomyocyte (CM) proliferation have been identified as potential candidates for functional heart repair following myocardial infarction. However, it is unclear whether the effects of these mitogens are species-specific and how they directly compare in the same cardiac setting. Here, we examined how CM-specific lentiviral expression of various candidate mitogens affects human induced pluripotent stem cell-derived CMs (hiPSC-CMs) and neonatal rat ventricular myocytes (NRVMs) in vitro. In 2D-cultured CMs from both species, and in highly mature 3D-engineered cardiac tissues generated from NRVMs, a constitutively-active mutant form of the human gene Erbb2 (cahErbb2) was the most potent tested mitogen. Persistent expression of cahErbb2 induced CM proliferation, sarcomere loss, and remodeling of tissue structure and function, which were attenuated by small molecule inhibitors of Erk signaling. These results suggest transient activation of Erbb2/Erk axis in cardiomyocytes as a potential strategy for regenerative heart repair.