1. Physics of Living Systems
  2. Structural Biology and Molecular Biophysics
Download icon

Dynamic metastable long-living droplets formed by sticker-spacer proteins

  1. Srivastav Ranganathan
  2. Eugene I Shakhnovich  Is a corresponding author
  1. Department of Chemistry and Chemical Biology, Harvard University, United States
Research Article
  • Cited 0
  • Views 529
  • Annotations
Cite this article as: eLife 2020;9:e56159 doi: 10.7554/eLife.56159


Multivalent biopolymers phase separate into membrane-less organelles (MLOs) which exhibit liquid-like behavior. Here, we explore formation of prototypical MOs from multivalent proteins on various time and length scales and show that the kinetically arrested metastable multi-droplet state is a dynamic outcome of the interplay between two competing processes: a diffusion-limited encounter between proteins, and the exhaustion of available valencies within smaller clusters. Clusters with satisfied valencies cannot coalesce readily, resulting in metastable, long-living droplets. In the regime of dense clusters akin to phase-separation, we observe co-existing assemblies, in contrast to the single, large equilibrium-like cluster. A system-spanning network encompassing all multivalent proteins was only observed at high concentrations and large interaction valencies. In the regime favoring large clusters, we observe a slow-down in the dynamics of the condensed phase, potentially resulting in loss of function. Therefore, metastability could be a hallmark of dynamic functional droplets formed by sticker-spacer proteins.


Biomolecular phase transitions are widespread in living systems. Transitions that result in irreversible, solid-like assemblies such as amyloid fibrils are a hallmark of disease while those like cytoskeletal filaments play a functional role (Banani et al., 2017). The third self-assembled state, in addition to the soluble and the solid-like state, is the loosely associated droplet phase held together by several weak, transient interactions (Guo and Shorter, 2015). The transient nature of these interactions makes these self-assemblies reversible and thereby a potential strategy for a temporally regulated sub-cellular organization. Several examples of spatiotemporally regulated, droplet-like objects within the cell have been discovered in the past few decades, composed of different types of proteins often co-localized with nucleic acids (Hyman et al., 2014). The importance of these membrane-less compartments is potentially two-fold: (i) localizing biochemical reactions within the cell, and (ii) sequestering biomolecules to regulate their activity (Hyman et al., 2014). Examples of cytoplasmic membrane-less organelles include P bodies, germ granules and stress granules (SGs) (Mitrea and Kriwacki, 2016). However, aberrant granule dynamics and a transition from a liquid-like to a more solid-like state are often hallmarks of degenerative diseases. (Kim et al., 2013; Ramaswami et al., 2013) Solid-like RNP aggregates have been reported as cytoplasmic inclusions (Patel et al., 2015) and as nuclear RNP aggregates (Ramaswami et al., 2013) in degenerative diseases. Investigating the physical principles governing the formation of these high-density phases is vital to understand the subcellular organization and the conditions leading to disease.

Several experimental studies have highlighted the ’multivalent’ nature of the constituent proteins in membrane-less organelles (Banani et al., 2017; Li et al., 2012; Xing et al., 2018). In other words, these proteins carry multiple associative or ’adhesive’ domains (Patel et al., 2015; Li et al., 2012). Multivalency could be achieved by several different architectures (Banani et al., 2017), the simplest being a linear sequence of folded domains that are connected by linker regions. Li et al., employed a linear multivalent 2-component model system (SH3 and PRM domains threaded together by flexible regions) to show that liquid-like droplet formation can result from just two multivalent interacting components (repeats of the same domain) (Li et al., 2012). Another intriguing feature of condensate proteins is the presence of intrinsically disordered regions or low complexity sequences that link folded domains together (Harmon et al., 2017). Therefore, a combination of these motifs and different architectures could form a basis for different types of phase-separated structures within the cell. In this paper, we address a specific question – how does this multi-domain sticker-spacer architecture shape the dynamics of droplet formation and further growth? Previous computational studies, notably the coarse-grained simulations by Harmon et al., address the nature of the equilibrium gel-like state formed by multivalent polymers with a sticker-spacer architecture (Harmon et al., 2017; Choi et al., 2019). The equilibrium phase diagrams for the 2-component SH3-PRM system and the effect of regulator molecules have also been studied using patchy particle models by Zhou and co-workers (Ghosh et al., 2019; Nguemaha and Zhou, 2018). Theoretical studies have also employed the Flory theory of phase separation in polymer solutions to describe the thermodynamics of membraneless organelle (Brangwynne et al., 2015) formation. The Flory theory, which applies to solutions of homopolymers predicts two phases – fully mixed and a single, large phase-separated droplet. Equilibrium theories robustly establish the underlying two-state equilibrium landscape that drives biopolymers to localize into ‘polymer-rich’ phases within the cell (Harmon et al., 2017; Choi et al., 2019; Brangwynne et al., 2015; Flory, 1942). Indeed, the growth of liquid droplets via the mechanisms of Ostwald ripening, coalescence or by consumption of free monomers from the surrounding medium (at concentrations greater than the saturation threshold) has been observed in vitro (Guillén-Boixet et al., 2020; Zhang et al., 2015; Martin et al., 2020). However, in a deviation from the equilibrium picture, membraneless organelles such as the nucleoli and P bodies are also known to exist in the form of multiple co-existing droplets (with equilibrium-like droplet environment) in vivo at biologically relevant time-scales (Brangwynne et al., 2011; Kilchert et al., 2010; Berciano et al., 2007). (Dine et al., 2018) report the persistence of multiple optogenically controlled droplets over long periods of time in vitro, a finding that is attributed to a progressive slow-down over time of droplet growth via Ostwald ripening. Active ATP-dependent processes have been attributed to the coexistence of multiple droplet phases, and the regulation of their size distributions in the cell (Wurtz and Lee, 2018a; Weber et al., 2019; Zwicker et al., 2017). Other potential factors hindering droplet fusion could be structural rigidity of the droplet scaffold (Boeynaems et al., 2019) and the aging of droplets due to internal re-organization (Lin et al., 2015). However, the potential existence of a more general, ATP-independent mechanism resulting in the prevalence of a stable multi-droplet system over biologically relevant time-scales (Wegmann et al., 2018; Li et al., 2012) has not been explored yet.

Despite significant experimental and computational efforts, a key question remains unanswered – what are the physical mechanisms that result in a long-living, metastable, multi-droplet system (Wegmann et al., 2018; Dine et al., 2018) in the absence of active processes? How does the sticker-spacer architecture of self-assembling biopolymers, with finite valencies, contribute to the slow-down and arrested coalescence of the early droplets into a single macro-phase? Since the equilibrium states in such systems are either mixed or a single, large droplet state, the kinetic factors must play a crucial role in the prevalence of the metastable multi-droplet state (Li et al., 2012; Dine et al., 2018). Also, if the system of multiple droplets is metastable, what are the physical factors that govern the kinetics of cluster growth? Measurable quantities such as droplet sizes also assume importance in the context of biological function, where the size of the membraneless organelles has been linked to function (Brangwynne, 2013). It is, therefore, not just the tendency to phase separate but also the size of these assemblies that must be regulated for their proper functioning (Goehring and Hyman, 2012). In this context, understanding the potential molecular mechanisms tuning the dynamics of cluster growth at biologically relevant time-scales becomes extremely relevant.

In this work, we address these questions using multi-scale coarse-grained models. In the first part of the paper, we probe the physical determinants of the formation of long-living metastable micro-droplets using a coarse-grained model of multivalent polymers composed of a linear chain of adhesive domains separated by semi-flexible linkers. The results of our Langevin dynamics (LD) simulations for this model shed light on the early stages of droplet growth and the mechanism of arrested phase-separation. Next, using the LD simulations as the basis to identify vital time-scales for condensate growth, we explore the phenomenon at biologically relevant time-scales using a phenomenological kinetic model. Broadly, we explore the kinetically arrested liquid-liquid phase separation (LLPS) on multiple scales – from mesoscale to macroscopic. Overall, these results provide a detailed mechanistic understanding of the factors that determine the prevalence of the metastable multi-droplet state for spacer-sticker heteropolymers in the absence of active, ATP-dependent processes. While the thermodynamic landscape drives phase-separation, the crucial role played by dynamic processes could result in measurable quantities deviating from the equilibrium predictions. The current study sheds light on the importance of dynamics of cluster growth, complementing the existing equilibrium understanding of the phenomenon.


Despite the complexity of the intracellular space, experiments suggest that in vitro liquid-liquid phase separation (LLPS) can be achieved even using simple two-component systems (Li et al., 2012). The tractability of simpler models makes them powerful tools to investigate the role of physical factors in modulating droplet formation and growth. Here, we perform LD simulations (see Materials and methods for detail) to understand the process of self-association between two types of polymer chains composed of specific interaction sites (red and yellow beads in Figure 1). These adhesive sites are linked together by non-specifically interacting linkers (blue beads in Figure 1). The red and yellow beads on these chains mimic complementary domains on different chains that can participate in a maximum of one specific interaction (between yellow and red beads).


Schematic of the polymer model for studying phase-separation by multivalent biopolymers.

For our simulations, we consider 400 such semi-flexible polymer chains (200 of each type) in a cubic box with periodic boundary conditions). Each chain in the simulation box is composed of 5 specific interaction sites that are linked together by non-specific linker regions that are 35 beads long. (blue beads in Figure 1). This linker length was based on previous theoretical studies of phase-separating proteins (Harmon et al., 2017). We employ conventional Langevin dynamics simulations to study the self-assembly of the model biopolymers, wherein the size of the specific interaction sites (diameter of 20 Å) is roughly four times that of the linker beads which represent individual amino acid residues (diameter of 4.2 Å). This difference in sizes is to mimic a folded adhesive domain – SH3 domain (diameter of 20 Å, PDB ID:1SHG,) that is often involved in liquid-liquid phase separation (Musacchio et al., 1992). The folded adhesive domains were modelled at a lower resolution (one bead per domain of 60 amino acids) than linkers which are more dynamic and hence modelled with a one bead per amino acid resolution.


Terminology and order parameters used in this study

Self-assembly of multivalent polymers with inter-chain interactions between complementary domains is a simple model system for studying phase-separation by intracellular polymers. In the current study, we employ this model to discover microscopic factors driving phase separation. In the first part of this study, we present results of self-assembly driven by irreversible (highly stable) functional interactions and identify two key time-scales that define cluster growth and their size distributions using Langevin dynamics (LD) simulations. In the latter part of this paper, we build a coarse-grained kinetic model demonstrating the tunability of this phenomenon using kinetic Monte Carlo (kMC) simulations. In the following subsections, we use the term specific interactions for the finite valency interactions between functional domains and non-specific interactions to refer to the inter-linker isotropic interactions. It must be noted that the emphasis of our simulations is to shed light on the factors that suppress the coalescence of individual clusters into a single, large equilibrium-like cluster. Also, the focus of the current study is to understand the process of cluster growth and the physical mechanisms modulating measurable quantities such as droplet size distributions at the early (LD) and biologically relevant time-scales (kMC). In Table 1, we summarize the key simulation variables and order parameters used in this study.

Table 1
Important simulation variables and order parameters.
Notation/TerminologyPhysical InterpretationDefinition
NtotTotal number of polymer chains in the system. Ntot = 400 in our simulations.
LclusSize of the single largest clusterrepresented as fraction of Ntot
SclusSize of the clusterrepresented as fraction of Ntot
CmonoConcentration of polymer chains in the simulation boxIn units of μM
φBulk density of proteins (in their monomeric state) when the individual chains are randomly placed in the simulation box at the start of the simulation.ϕ=Ntot(4/3)πRG3,.Rg is the the radius of gyration of proteins when they are
randomly positioned in simulation box at the start of the simulation.
φclusIntra-cluster density of polymer chains.ϕclus=Sclus(4/3)πRg3clus.Rgclus is the radius of gyration of the system of proteins within the cluster.
φclusNormalized intracluster density describing the degree of enrichment of polymer chains within the cluster.For system-spanning networks,
φclus /φ→1. For dense clusters,
φclus /φ >> 1.
εnsInteraction strength for isotropic, non-specific interactions between linker regions. εns in LD simulations is a pairwise interaction strength between individual beads In kMC simulations, εns is the net non-specific interaction strength between two lattice particles.

εspStrength of attractive interaction between functional domains (specific interactions).
φlatticeBulk density of monomers in the 2D-lattice in kMC simulations. This quantity is analogous to the concentration of monomers on the lattice.φlattice = Ntot/L2, where L is the size of the 2D square-lattice.
λValency of the polymer chain, that is the number of adhesive functional domains per interacting polymer chain.
ΚBending stiffness of the linker regions in LD simulations. A higher value of K is used to model stiffer linkers that prefer a more open configuration.Κ = 2 kcal/mol, in LD simulations with flexible linkers.
ηThe viscosity of the medium in LD simulationsη = 10−3Pa.s in LD simulations, unless mentioned otherwise.
kdiffDiffusion rate of free monomers in the 2D-lattice kMC simulations.
kbondRate of formation of specific interactions between neighboring particles in 2D-lattice kMC simulations.

When LclusNtot—> Ntot, we refer to the state as a system-spanning ‘macro-phase’. A system of multiple, co-existing clusters with with SclusNtot<< Ntot, is referred to as a metastable ‘micro-phase’. It is important to note that the use of the term ‘micro-phase’ in this paper is in reference to metastable, and not equilibrium phases. In order to establish the condensed nature of micro-phase clusters as compared to the system-spanning macro-phase, we characterize the intra-cluster densities normalized by the density of a system (ϕclus /ϕ) of randomly placed polymer chains in the simulation box. We use an analogous term for density measurements in our kinetic Monte Carlo simulations, with the volume term being replaced by the area of the cluster (or lattice) in 2D. For details of the potential functions and the implementation of these simulations, please refer to the Materials and methods section.

Metastable ‘micro-phases’ – a likely outcome at low to intermediate concentrations

Multivalent polymers with adhesive domains separated by flexible linkers can potentially self-assemble through two types of interactions, a) the finite number of adhesive contacts between the functional domains (yellow and red beads in Figure 1), and b)non-specific, isotropic interactions between the linker regions (blue beads in Figure 1). We first performed control simulations with specific interactions turned off, where we varied free monomer concentration Cmono, from 10 to 200 μM for a weak non-specific interaction strength of ϵns = 0.1 kcal/mol. In these control simulations, we observe no phase-separation in the whole range of Cmono(Figure 2—figure supplement 1A, purple curve) . Therefore, in this regime of ϵns and Cmono, at the simulation time-scale of 16 μs, the polymer assemblies do not reach large sizes. Not surprisingly, this result suggests that the assembly driven by non-specific interactions (linker-driven) alone is achieved only at stronger non-specific interactions and/or high free monomer concentrations. However, intracellular phase-separation often results in the enrichment of biomolecules within the self-assembled phase, at relatively low protein concentrations (Xing et al., 2018). Under such conditions, specific interactions which are fewer in number become critical determinants of phase separation. In our LD simulations, we employ polymer chains comprised of 5 univalent adhesive domains and four linker regions totalling 145 beads, bringing a total valency of 5 for each polymer. For the sake of simplicity, we begin with a situation where these bonds, once formed, do not break for the rest of the simulation. Although an idealized construct with respect to biological polymers, these simulations can build intuition about cluster size distributions and arrest of cluster growth when the spontaneously assembled clusters are unable to undergo further re-organization.

In Figure 2A we show the single largest cluster size, for varying free monomer concentrations (Cmono). The polymer-chains, for the data plotted in Figure 2, have flexible linkers (κ=2 kcal/mol, see Table 1. and Materials and methods for definition of bending rigidity) with weak inter-linker interactions (ϵns=0.1 kcal/mol, see Table 1. And Materials and methods section for definition). As evident from Figure 2, the polymer chains can assemble into larger clusters in the presence of functional interactions (black curve, Figure 2A), as opposed to the control simulations where inter-linker interactions alone drive self-assembly (Figure 2—figure supplement 1A, purple curve). Simulations with irreversible functional interactions suggest that, for low and intermediate concentrations (Cmono<70 μM), SclusNtot<< Ntot. Also, at these concentrations, we observe a distribution of cluster sizes (Figure 2C, blue bars) indicating a multi-cluster ‘microphase’ state. We also characterized the density of the clusters using a second order parameter – ϕclus /ϕ. This quantity is analogous to the measure of phase separation employed in a previous study by Harmon et al., 2017. The density transition as a function of concentration shows a non-monotonic behavior. For large Cmono, where the mean largest cluster sizes approach the size of the system, we do not observe a polymer-dense phase. In this regime (pink region in Figure 2A, where ϕclus/ϕ1, the self-assembled state is a percolated, system spanning network. In a narrow intermediate regime of concentration (cyan region in Figure 2A), we observe microphases with high polymer density within the cluster. The density transitions are qualitatively consistent with predictions from equilibrium theories for sticker-spacer proteins studied by Harmon et al., 2017. Although the density transitions are in tune with the equilibrium landscape, the key prediction of our simulations is the microphasic nature of this high-density regime with coexistence of multiple clusters (Figure 2A and C). The intra-cluster density, however, would also depend on the characteristics of the linker, a phenomenon we demonstrate in subsequent sections.

Figure 2 with 1 supplement see all
Cluster Sizes in Langevin Dynamics simulations.

(A) Black Curve – The single largest cluster as a function of free monomer concentration (in μM) for irreversible functional interactions. The largest cluster size is shown as a fraction of the total number of monomers in the simulation box. Purple Curve – The mean density of the cluster, ϕclus = Sclus(4/3)πRg3clus, normalized by the bulk density of monomers in the simulation, ϕ = Sclus(4/3)πRg3. Sclus and Ntot refer to the size of the cluster and the total number of monomers in the simulation box, respectively. Rgclus and Rg refer to the radius of gyration of the cluster, and the radius of gyration of proteins in their randomly located initial configuration at the start of the simulation, respectively. The quantity ϕclus/ϕ shows the degree of enrichment of polymer chains within the cluster upon self-assembly. The smooth curves are plotted as a guide to the eye, using the cspline curve fitting. (B) Comparison of the mean sizes of the single largest cluster for reversible and irreversible specific interactions, for varying free monomer concentrations. (C) Cluster size distributions for varying free monomer concentrations for reversible and irreversible specific interactions. Darker shades of red in the distributions indicate verlapping regions of the distribution while the blue and light red shades indicate regions with no overlap in the presence and absence of breakable interactions. The linker stiffness for the self-assembling polymer chains in this plot is 2 kcal/mol while the strength of inter-linker interaction is 0.1 kcal/mol (per pair of interacting beads). The mean and distributions of the largest cluster sizes were computed using 500 different configurations from five independent simulation runs of 16 μs.

To further test whether this early time-scale behavior changes in the presence of reversible interactions, we introduced breakable specific bonds in our model (see Materials and methods section for details of the implementation). As with the irreversible interaction simulations, even in the presence of breakable interactions, LclusNtot << Ntot (Figure 2B) except for large Cmono when we observe a system-spanning network. Critically, we find coexistence of intermediate cluster sizes with small and large clusters suggesting an increased diversity in cluster sizes at the early stages of droplet assembly upon introduction of breakable interactions (Figure 2C).

Overall, these results suggest that in the concentration regime of dense clusters (Figure 2A, blue shaded region), we observe a distribution of cluster sizes (Figure 2C) as opposed to a single, large cluster. Given that the cellular concentration of phase-separating proteins is often in the nanomolar to the low micromolar range (Xing et al., 2018), the multi-droplet state could persist in vitro and in vivo even in the absence of active processes. In the subsequent sections, we explore the mechanisms that could prevent or significantly delay the coalescence of multiple small clusters into a single large equilibrium-like cluster, eventually resulting in a distribution of droplet sizes.

Exhaustion of free valencies results in kinetically arrested droplets

Our LD simulations so far reveal an interesting trend – while the intracluster density transition shows a non-monotonic dependence on concentration that is qualitatively consistent with equilibrium theoretical predictions (Harmon et al., 2017), we do not observe a single, large equilibrium-like cluster at low concentrations. Crucially, in the regime where the intra-cluster densities are at their highest, we observe a distribution of cluster sizes (Figure 2A and C), with SclusNtot << Ntot, suggesting the prevalence of arrested micro-phases. Therefore, identifying the time-scales which are vital for cluster growth could reveal the cause of arrested droplets growth. In Figure 3A and B, we show the time evolution of the individual aggregate species at different Cmono (10μM and 50 μM). As seen from Figure 3A and B, for irreversible functional interactions, the monomer fraction continues to monotonically decrease during the simulations with the fraction of other competing species increasing concomitantly. However, at low concentrations (Figure 3A, 10 μM), the monomer fraction curve (Figure 3A, purple curve) shows a cross-over with the dimer curve (Figure 3A, green curve) while higher-order clusters do not appear at simulation time-scales. This result suggests that the spontaneous formation of large assemblies held together by functional interactions is contingent upon two time-scales. The first is the diffusion-limited time-scale that governs initial dimerization and the subsequent growth of these smaller clusters. The second, competing time-scale is the one where all functional valencies get exhausted within the smaller initial clusters. At an intermediate concentration of 50μM, (Figure 3B) we observe that the fraction of large aggregates (>10 mer) increases during the early part of the simulations before free monomers are entirely consumed within dimers. In other words, the unsatisfied valencies within the monomers, dimers and trimers get utilized to form larger-sized clusters before the specific valencies get exhausted within the smaller clusters making them no longer available for further self-assembly. To further verify this observation, we track the temporal evolution of the fraction of available valencies within and outside the single largest cluster (Figure 3C). We observe that while there are unutilized valencies within the single largest cluster (Figure 3C, purple curve), the cluster does not grow further due to almost complete exhaustion of valencies within the polymer chains outside the single largest cluster (Figure 3C, green curve). A slower rate of exhaustion of free-valencies within micro-clusters by varying the specific bond formation probability, Pform , results in larger cluster sizes for smaller free monomer concentrations (Figure 3—figure supplement 1), suggesting that the dynamicity of bond formation can alter the cluster size distributions significantly. These results suggest that the ability to form a single, large, macro-cluster is limited by the exhaustion of free valencies within smaller sized clusters, thereby arresting their growth (Figure 3D).

Figure 3 with 1 supplement see all
Tracking cluster formation at early timescales.

A and B show the temporal evolution of specific contacts for a free monomer concentration of 10, and 50 μM, respectively. For a low concentration of 10 μM, there is an initial decrease in the monomer population (purple curve) which is concomitant with an increase in the dimer population. A negligible fraction of the clusters is in the form of large-mers (size > = 10, orange curve) at these low concentration since the available valencies for growth are consumed by the smaller aggregated species. An increase in concentration from 10μM to 50μM results in an increase in the large-mer (orange curve, 50 μM) population as the monomer fraction decreases during the simulation. A higher free monomer concentration allows the larger clusters to grow due to consumption of free monomers (with unsatisfied valencies) before they get converted into smaller clusters (dimers, trimers) with satisfied valencies. (C) Time evolution of available valencies within the single largest cluster and outside the single largest cluster, for a Cmono of 50 μM, ϵns of 0.1 kcal/mol and linker bending rigidity of 2 kcal/mol. (D) A schematic figure showing the possible mechanisms of cluster growth and arrest and the competing timescales that could punctuate the process.

Identifying the vital timescales determining cluster growth

For stable functional interactions, the process of phase separation gets arrested due to kinetically trapped clusters which do not participate in further cluster growth due to lack of available valencies (Figure 3D). As discussed above, two critical time-scales would then dictate the growth of clusters: i) the time-scale for two chains to meet and form the first functional interaction, and ii) the time it takes for the polymer chains within an assembly to exhaust all valencies within the cluster before new chains join in. We explore the factors that these two timescales depend on, using the primary unit of any self-assembly process – the dimer. Using LD simulations, we first compute the mean first passage times for two polymer chains to form their first functional interaction. Given the diffusion-limited nature, this time-scale varies inversely with the concentration of free monomers (represented in the form of density) in the system (Figure 4A).

Factors influencing the key timescales.

(A) The mean first passage time for the first specific interaction between a pair of polymer chains as a function of the bulk density, ϕ, (see Table 1 for definition). (B) The mean first passage time for a pair of polymer chains to exhaust all the available valencies within the dimer.

This diffusion-limited time-scale dictates the encounter probability of the two chains. Once the first bond is formed, resulting in an ’active dimer’ (one that still has unsatisfied valencies), the cluster can only grow to larger sizes as long as the dimer remains active. Therefore, the time taken by the dimer to exhaust all its valencies becomes a vital second time-scale. In Figure 4B, we plot the mean first passage times for a dimer to exhaust all its valencies once the first bond is formed. For low linker stiffness (κ<2 kcal/mol in Figure 4B), the linkers behave like flexible polymers. In this limit (κ<2 kcal/mol), the mean time for the valency to exhaust within the dimer (referred to as the re-organization time, Texhaust-valency ) is independent of linker stiffness. For (κ>2 kcal/mol), the linkers behave like semi-flexible polymers and the re-organization time shows a dependence on the stiffness of the linker, with an increase in Texhaust-valency with κ. As the linker region becomes more rigid, this time-scale becomes slower, resulting in the dimer remaining ‘active’ for much longer. A slower re-organization time and a faster diffusion encounter time favors the cluster growing into larger sizes. It must however be noted that an increase in Texhaust-valency via increased linker rigidity would not only influence the cluster sizes but also the intracluster density. For stiffer linkers (and also more open linker configurations), the self-assembled state would approach a system-spanning network as shown in previous studies with linkers with greater effective solvation volumes (Harmon et al., 2017). Conversely, a faster re-organization time that is of the order of the encounter time-scales results in a higher likelihood of the clusters getting locked at smaller sizes. For stable functional interactions, these two time-scales would dictate the size of the largest cluster.

Linker regions as modulators of self-assembly propensity

The spatial separation of specifically interacting domains (with finite valencies) and the non-specific linker regions is an architecture that is highly amenable to being tuned for phase-separation propensities. In this context, we further extend the findings from our LD simulations to probe how the microscopic properties of the linker region could modulate the extent of phase-separation. In Figure 5, we demonstrate how the linker properties could be useful modulators of cluster sizes, without any alteration of the nature of the specific functional interactions (for a Cmono of 50μM). Further, to model an unstructured linker, we consider a scenario where the linkers participate in inter-linker interactions alone. Strong inter-linker interactions increased occupied valencies (Figure 5—figure supplement 1), for all the concentrations under study, while resulting in much smaller mean largest cluster sizes (Figure 5A) during the time-scales accessed by our simulations. Further, the mean radius of gyration (<Rgchain >) for the monomers within these clusters shows a sharp transition with an increase in inter-linker interaction strength (ϵns ). An increase in ϵns from 0.3 kcal/mol to 0.5 kcal/mol results in a sharp decrease in <Rgchain >, indicating an onset of linker-driven coil-globule transitions for the polymers within clusters (Lifshitz et al., 1978). This effect is also manifest in the increase in the density of the self-assembled clusters upon an increase in ϵns. For a Cmono of 50μM (Figure 5B), we see an increase in cluster density in the range of ≈5 to 50 times of that of the bulk upon variation in linker interactions. A similar trend was observed for a lower Cmono of 30μM, with a 10–100 fold variation in degree of enrichment within the cluster upon varying inter-linker interactions (Figure 5—figure supplement 2). Therefore, the density transitions could be tuned by varying the free monomer concentrations as well as the properties of the linker. This is consistent with previous equilibrium predictions (Harmon et al., 2017; Choi et al., 2019) that establish the dependence of density transitions on the properties of intrinsically disordered linkers. The mean density of clusters is also consistent with experimental findings of a ≈10–100 fold enrichment of biomolecules within droplets (Mitrea and Kriwacki, 2016; Li et al., 2012; Xing et al., 2018; Nott et al., 2015; Burke et al., 2015).

Figure 5 with 3 supplements see all
Linkers as modulators of self-assembly propensity.

(A) The size of the largest cluster for flexible linker regions(κ=2 kcal/mol) with varying inter-linker interaction strength (black curve, 0.1 kcal/mol and purple curve, 0.5 kcal/mol). Sticky inter-linker interactions result in smaller cluster sizes. (B). Purple Curve – Mean Radius of gyration for the individual-polymer chains (<Rgchain>) within a self-assembled cluster as a function of increased inter-linker interactions. Black Curve— Mean density of the clusters as a function of inter-linker interaction strength. (C) The mean re-organization times, Texhaust-valency as a function of linker-stiffness, for different values of inter-linker interaction strengths.

The assemblies that ensue at stronger inter-linker attraction are more compact and condensed akin to homopolymer globules (Lifshitz et al., 1978). We further probe the manner in which linkers can tune the dynamics of self-assembly, in addition to influencing the equilibrium behavior such as intracluster density. Intramolecular compaction of polymers due to non-specific inter-linker interactions brings specific domains closer in space, leading to a higher likelihood of the exhaustion of specific interaction valencies within small assemblies. In Figure 5C, we show how the inter-linker interaction strength can influence the time it takes to exhaust specific interaction valencies within a dimeric cluster (Texhaust-valency ). For weaker inter-linker attraction (ϵns<0.5 kcal/mol), the initial polymer assemblies are less compact (Figure 5B) and thereby exhaust valencies within a cluster at a much slower rate (higher Texhaust-valency for dimers with ϵns < 0.5 kcal/mol in Figure 5C). An increase in inter-linker interaction propensity results in faster re-organization times for these polymers. The polymers with ϵns = 0.5 kcal/mol exhaust their valencies almost an order of magnitude faster than their 0.1 kcal/mol counterparts (Figure 5C). Upon exhaustion of these specific interaction valencies, these clusters can only grow via inter-linker interactions, a phenomenon that could be less dynamic and tunable than the functional interaction driven cluster growth. It must, however, be noted that the observation of further coalescence of clusters formed by sticky inter-linker interactions was limited by the time-scales accessible to the LD simulations. Any alteration to the ’stickiness’ of the linker can shift the mechanism of assembly, and thereby result in altered kinetics of cluster growth by modulating the Texhaust-valency. Our dimerization simulations show that a second mechanism of slowing down the Texhaust-valency time-scale is by altering the flexibility of the linker region (increasing linker stiffness in Figure 5C). Primarily the stiffer linkers lead to more ‘open’ configurations of spacer regions. This is corroborated by a shift in the cluster size distribution towards larger sizes, for the polymer chains with rigid linkers (Figure 5—figure supplement 3). The linker region can thereby serve as a modulator of phase-separation propensity (characterized by the density) and cluster sizes. A variation in the intrinsic properties of the linker, with no modifications to the functional region, can be used as a handle to tune the density of polymers within the condensed phase (Figure 2 and Figure 5). This lends modularity and functionality to these condensates, with the linker regions influencing both the degree of enrichment and the dynamics of self-assembly.

Coarse-grained kinetic simulations predict the prevalence of metastable micro-phases at biologically relevant time-scales

The LD simulations help us identify the initial events that mark phase-separation by multivalent polymer chains assembling via finite-valency, specific interactions. However, the model is limited in its ability to access longer, biologically relevant time-scales at which droplets typically form and grow in living cells. We further explore the coalescence of multi-cluster system into a single, large equilibrium-like cluster at longer time-scales that are inaccessible to LD simulations. Here, we employ a coarse-grained approach wherein the whole polymer chain from the bead-spring model (with a fixed valency) is represented as an individual particle on a 2D-lattice. In our simulations, the lattice is populated by Ntot such multivalent particles (at varying densities, ϕlattice ) that diffuse freely, at a rate kdiff, and the number of particle collisions per unit time is proportional to kdiff*ϕlattice . In our kMC study, we vary ϕlattice in a range of 0.01 to 0.1 (see Materials and methods section for rationale). When two such particles occupy neighboring sites on the lattice, they interact non-specifically with an interaction strength of ϵns , a parameter that is analogous to the inter-linker interactions in our LD simulations (see Table 1). Additionally, two neighboring particles with unfulfilled valencies can form a specific bond (with finite valencies per lattice particle) with a rate of kbond. The valency per particle (λ) here can be utilized to form bonds with one to four potential neighbors, with each pair being part of one, or more than one specific interactions between themselves. However, unlike the irreversible specific interactions in our LD simulations, the specific bonds in the lattice model can break at a rate kbreak = kbond *exp(-ϵsp), where ϵsp is the strength of each specific bond. Additionally, clusters can diffuse with a scaled diffusion rate that is inversely proportional to the cluster size (kdiff /Sclus). It must be noted that the time-scale for the first bond formation in the LD simulations is an outcome of two phenomenological rates in this model, kdiff and kbond. The second time-scale, Texhaust-valency, is a time-scale that depends on the kbond and kbreak parameters in this model. The details of the simulation technique and the various rate processes can be found in the Materials and methods section and described schematically in Figure 6.

A schematic figure detailing the different rates in our phenomenological kinetic model simulated using the Gillespie algorithm.

The particles on the lattice can diffuse freely (when there are no neighboring particles) with a rate kdiff. In the presence of a neighboring particle, a non-specifically interacting monomer can diffuse away with a rate kdiffϵns.exp(- ϵns). Neighboring particles can also form specific interactions (with fixed valency λ) at a rate kbond or break an existing interaction with a rate kbreak. Clusters could diffuse at a rate that is scaled by their sizes Sclus. ϵns and ϵsp refer to the strength of non-specific and specific interactions, respectively.

Using kinetic Monte Carlo simulations (see Materials and methods and Gillespie, 1977), we explore the cluster formation (at times reaching a physiologically relevant scale of hours) by varying parameters such as, (a) bulk density of particles on the lattice (ϕlattice), (b) rate of bond formation (kbond), (c) valency per interacting particle (λ), and (d) the strength of specific interactions (ϵsp). With the assumption that, for diffusion-limited self-assembly, the rate of free diffusion kdiff is the fundamental time-scale limiting cluster growth, we first explore the relationship between the rate of specific bond formation (kbond), and kdiff (Figure 7A and Figure 7—figure supplement 1). It must be noted that the LD simulations employed the assumption that bond formation, upon the two functional domains coming in contact, is an instantaneous event. Here, we show that for values of kbond/kdiff0, there is no phase-separation. As the bond formation rate approaches that of free diffusion – corresponding to the instant bond formation assumption in LD simulations – we encountered phase-separated states in our simulations. However, the system largely favors the micro-phase separated state (bluish-red regions in the phase diagram, Lclus<<1) for low-intermediate ϕlattice even at biologically relevant time-scales of hours in the kMC simulations. In this regime of low-intermediate ϕlattice (bluish-red regions of the kinetic phase diagrams in Figure 7), we see denser clusters (Figure 7—figure supplement 2) with a wide cluster size distribution (distributions labeled μ1 in Figure 7—figure supplement 3). At high values of ϕlattice , we observe a system-spanning macrophase with lower cluster density (Figure 7—figure supplement 2) and a cluster size distribution that favors extremely large clusters co-existing with very small clusters (distributions labeled μ2 in Figure 7—figure supplement 3).

Figure 7 with 6 supplements see all
Kinetic Monte Carlo Simulations.

(A) Phase diagram highlighting the different phases (metastable microphase (μ1) or system-spanning macrophase (μ2), and the non-phase separated state (No PS) ) encountered upon increasing kbond(between 0 and kdiff ), and the bulk density of monomers (ϕlattice ) within the box. The assembling particles have a valency (λ) of 5 in these simulations. The region shaded yellow represents part of the phase-space where we observe a system-spanning network whereas the bluish-red region represents the metastable micro-phase with a distribution of cluster sizes (see Figure 7—figure supplement 3 for distributions). As we move from the blue to the red regions in this dynamic phase diagram, the size of the single largest cluster increases (also see Figure 7—figure supplement 1). (B) Phase diagram highlighting the different phases encountered upon varying valency and bulk density as the phase parameters. The system-spanning macrophase is only encountered for larger valency particles at higher densities. This phase diagram was computed for kbond=kdiff and ϵns of 0.35 kT. (C) The fraction of monomers in the largest cluster as a function of epsilon, for kdiff = kbond and ϕlattice = 0.04. The single largest cluster sizes in all sub-panels of this Figure were computed for a simulation time-scale of 2 hr (with the fundamental timsecale of diffusion being set to kdiff = 1 s-1 ).

Further, this phase diagram (Figure 7A) also establishes that for the value of non-specific interaction strength (ϵns=0.35 kT, see Materials and methods section for rationale) used here the cluster formation is driven by the finite-valency specific interactions (absence of clustering for kbond/kdiff0). For comparison, in Figure 7—figure supplement 4, we present the mean cluster sizes for assemblies that are stabilized by non-specific interactions only (λ=0). The ϵns-ϕlattice phase diagram shows that non-specific interaction-driven cluster formation occurs at only high values of ϵns (Figure 7—figure supplement 4B) . Therefore, cluster formation is contingent on the bonding rate being of the same order as the free-diffusion rate (kbondkdiff) establishing the validity of the instantaneous bonding assumption in the LD simulations. It is the ratio of kbond/kdiff , and not the absolute magnitudes, that is a vital parameter for these simulations. Hence, in all the kMC simulations we set the value of kdiff to 1 s-1 and vary the ratio kbond/kdiff to tune phase separation. It must be noted that, unless mentioned otherwise, the results from the kMC simulations presented here are for a weak non-specific interaction strength of 0.35 kT. All simulations were performed for a time-scale of 2 hr (actual time). As proof of convergence of these simulations, we compare results at the end of 2 hr to those at longer simulation time-scale and show that there is negligible difference in cluster sizes (Figure 7—figure supplement 5).

We systematically explored the effect of valency of specific interactions on the extent of phase separation. Figure 7B shows kinetic phase diagram with λ and ϕlattice as the phase parameters. For smaller λ and low ϕlattice , Lclus (and Sclus ) NtotNtot (blue and black regions in Figure 7B, and Figure 7—figure supplements 1 and 3). This suggests that, at biologically relevant concentrations and time-scales, exhaustion of free valencies for particles with smaller λ is a crucial determinant of cluster sizes. This phase-diagram is consistent with in vitro experiments involving SH3 and PRM chains with varying valencies, with higher valency molecules displaying a lower critical concentration for phase separation (Li et al., 2012).

In addition to the valency of particles, a vital parameter that would determine the droplet sizes is the strength of these specific interactions (ϵsp ). As evident from Figure 7C, for specific interactions that are extremely weak (ϵsp<2 kT in Figure 7C), there is no significant phase separation. Strikingly, this critical interaction strength (when the largest cluster >= 50% of available monomers) is lower for higher valency particles (λ=5 curve in Figure 7C). Interestingly, the SH3-PRM interaction strength is reported to be in the range of 2–5 kT (Li et al., 2012; Harmon et al., 2017). A more detailed ϵsp -kbond kinetic phase diagram can be found in Figure 7—figure supplement 6. The properties of the condensate at biologically relevant time-scales could thus be tuned via different parameters, offering the cell several handles to modulate sizes and morphologies of droplets.

Tunability of exchange times for the metastable micro-droplets

In the kMC simulations so far, we discuss the manner in which different phase parameters could shape droplet size distributions. However, the functionality of a condensate hinges not only on the ability of biomolecules to assemble into larger clusters but also to exchange components with the surrounding medium at biologically relevant time-scales. These exchange time-scales are also a measure of the material properties of the droplets themselves (Guo and Shorter, 2015; Xing et al., 2018; Feric et al., 2016). Therefore a systematic understanding of the dependence of molecular exchange times on intrinsic and extrinsic parameters is crucial to get a grasp of the tunability of intracellular self-organization driven by finite-valency specific interactions. In this context, we probed the extent to which the exchange times could be tuned by modulating the intrinsic features of the self-assembling units. Here, we define monomer exchange times as the mean first passage time for a monomer to go from having four neighbors to being completely free. To compute first passage times, we kept track of exchange events from across 100 simulation trajectories of 10 hr each. Our simulations suggest that, for a given valency, a slight increase in interaction strength ϵsp within a narrow window could result in a dramatic increase in the size of the clusters (Figure 7D). This raises an interesting question—is there an optimal range for these phase parameters that promote phase separation while maintaining the dynamicity of the clusters? In this context, we first computed the mean first passage times for monomer exchange upon systematic variation of ϵsp (Figure 8A). Interestingly, as with cluster sizes, a slight increase in ϵsp results in a dramatic increase in monomer exchange times. For particles with λ=5, a slight increase in ϵsp from 2 kT to 2.5 kT, there is a four-fold slow-down in exchange times indicating dramatic malleability in the dynamicity of these assemblies. This shift from the fast to slow exchange dynamics is less abrupt in case of particles with λ=3 suggesting that an interplay between λ and ϵsp could tune the droplets for desired exchange properties. We further varied these two parameters (λ and ϵsp) systematically to probe their effect on cluster sizes (Figure 8—figure supplement 1A) and molecular exchange times (Figure 8—figure supplement 1B). As expected, for weak ϵsp and a low λ there is no phase separation (black region in Figure 8B). For an intermediate regime in this phase-space, the system is predominantly in a metastable micro-phase separated state, with either slow (red region in Figure 8B) or fast (blue region in Figure 8B) molecular-exchange times. However, a system-spanning macro-state is only observed for really large λ and ϵsp, with a dramatic slow-down in exchange times (shaded yellow region in Figure 8B and Figure 8—figure supplement 1B), suggesting that these assemblies might be biologically non-functional. Valency is, therefore, a key determinant of how frequently a molecule gets exchanged between the droplet and the free medium. Similar observations have been made experimentally by Xing et al. (Figure 8C and D) with regards to several condensate proteins featuring different valencies, with low valency species showing exchange times that are orders of magnitude faster than the higher valency ones (Xing et al., 2018). Given that functional droplets are tuned for liquid-like behaviour, metastable microphases are, therefore, a likely outcome for dynamically exchanging droplets.

Figure 8 with 1 supplement see all
Effect of model parameters on the exchange times between monomers and the aggregates.

The parameter values used in panels A and B are ϕlattice=0.04, kbond/kdiff = 1. (A) Mean first passage time for the monomers to go from the buried state (with four neighbors) to the free state (with no neighbors) in response to varying values of specific interaction strength. The two curves show the trends for species with different interaction valencies. The dashed solid lines refer to value of ϵsp beyond which the largest cluster consumes 50% or more of available monomers. (B) The state of the system, for variation in ϵsp and λ suggests that the system displays remarkable malleability in dynamicity and size distributions. Weak ϵsp and a low λ results in no phase separation (shaded black). For higher values of both these parameters, the system can access the single large macro-phasic state (shaded yellow), however with a dramatic slowdown in exchange times. For an intermediate range, the system adopts a metastable microphase separated state, with either slow (shaded red) or fast-exchange (shaded blue) dynamics. (C) The experimentally determined molecular exchange times for molecules of varying interaction valencies. (Xing et al., 2018). (D) The extent of recovery after a photobeaching experiment, for interacting species with varying valencies. The data in panel C) and D) has been obtained from a study by Xing et al., 2018. The solid curves in C and D are added to guide the eye.

Effect of solvent viscosity on dynamics of cluster growth

The coarse-grained LD and kMC simulations suggest that kinetic barriers to cluster growth result in the frequently observed state of the long-living multi-droplet state in vitro and in vivo. To highlight this further, we performed LD and kMC simulations where we vary the diffusive properties of the free monomers in solution by tuning solvent viscosity (LD) and diffusion rate (kMC simulations). We first performed LD simulations for solvent viscosities that are 1.5 and 3 times that of water (η=10-3 Pa .s). As evident from the cluster size distributions at the end of a 20 μs simulation, the larger clusters were not observed for more viscous solvents (shaded region in Figure 9A and B). This is further corroborated by the smaller cluster sizes at higher solvent viscosities (Figure 9C). This suggests that slowing down the diffusion-limited time-scale for higher viscosities could result in a further reduction in the observable cluster sizes at finite time-scales relevant to biology. At higher viscosities, we observe large clusters only at higher free monomer concentrations (Figure 9—figure supplement 1). We further probed the role of diffusion-limited time-scale on the observed cluster sizes at a fixed time-scale of 2 hr in our kinetic Monte Carlo simulations. As with our LD simulations, slower diffusion rates result in smaller cluster sizes for a simulation time-scale of 2 hr, both for a λ of 3 and 5 (Figure 9D and E). These results are consistent with previous in vitro studies which suggest that solvent viscosity could influence the rate of droplet coalescence (Kaur et al., 2019). These results make a key testable prediction of this work pointing to the vital role of dynamics in the formation of membraneless organelles as metastable micro-droplets positing that the measurable quantities such as droplet sizes depend on dynamic quantities of the solvent such as viscosity. If micro-droplets were formed at equilibrium, such dependencies would not be observed (Harmon et al., 2017; Brangwynne et al., 2015). While initial indications point out to the crucial role of solvent viscosity (Kaur et al., 2019) more experimental studies are needed to confirm or refute this key prediction of our study.

Figure 9 with 1 supplement see all
Effect of solvent viscosity.

Langevin Dynamics Simulations. (A and B) Cluster size distributions from Langevin dynamics simulations for different values of solvent viscosity and different free monomer concentrations. The shaded region in the two subpanels (A) and (B) is to highlight that a more viscous solvent results in larger clusters disappearing from the distribution. (C) The size of the single largest cluster as a function of solvent viscosity for a free monomer concentration of 70 μM. (D and E) Monte Carlo Simulations. The size of the single largest cluster as a function of density for varying rates of diffusion, for specific interaction valencies λ=3 (D) and λ = 5 (E).


Long-living metastable droplets: a potential signature of multivalent heteropolymers

Membrane-less organelles are heterogeneous pools of biomolecules which localize a high density of proteins and nucleic acids (Banani et al., 2017; Feric et al., 2016; Boeynaems et al., 2019). An interesting feature of several complex macromolecules (Mitrea and Kriwacki, 2016; Li et al., 2012; Harmon et al., 2017) that constitute these droplets is ’multivalency’ stemming from multiple repeats of adhesive domains (Banani et al., 2017; Li et al., 2012; Patel et al., 2015). These adhesive domains can bind to complementary domains on other chains, thereby facilitating phase separation. In this study, we model this phenomenon as that of self-associative polymers that possess folded domains (represented as idealized spheres in Figure 1) separated by flexible linker regions. Recent computational studies have employed similar models to characterize the equilibrium state of these polymer systems, notably the coarse-grained simulations by Harmon et al., 2017 and Choi et al., 2019. These studies employ the ’sticker and spacer’ model to understand the phase behaviour of linear multivalent polymers (Harmon et al., 2017), mainly focusing on the nature of the phase-separated state at equilibrium. Using lattice-polymer Monte Carlo simulations, they establish the role of intrinsically disordered linker regions as molecular determinants that dictate the equilibrium state of a system of associative polymers that interact via non-covalent interactions (Harmon et al., 2017). Crucially, these works focus on the cross-linked gel-like nature of the equilibrium state of these polymers and establish the underlying thermodynamic landscape governing phase separation of spacer-sticker polymers. The observation of the single, large equilibrium-like droplet phase in vitro establishes the robustness of the existing thermodynamic understanding of the phenomenon. However, experiments also report the presence of multiple, co-existing droplets (with equilibrium-like droplet micro-environments) that remain stable at biologically relevant time-scales (Brangwynne et al., 2011; Kilchert et al., 2010; Dine et al., 2018; Wegmann et al., 2018) in vivo and in vitro. Minimal models by Wurtz and Lee, 2018b point towards the potential role of ATP-dependent processes (Weber et al., 2019) in suppressing Ostwald ripening, and thereby the stabilization of the multi-droplet system. However, the potential role of a more general, ATP-independent mechanism in establishing a long-living multi-droplet phase has not yet been explored.

The deviation from the equilibrium picture raises a number of interesting questions –how does the multivalent, multi-domain architecture of proteins influence the suppression of droplet growth into a single macroscopic phase by coalescence? Does the physics of associative polymers such as multivalent proteins differ from those of simple, homopolymer chains? Also, in the context of the size-dependent function of membrane-less organelles, can the intrinsic features of the self-assembling chain act as handles that control the dynamics and size distributions of the droplets? In this paper, we argue that, even in the absence of active processes, the finite interaction valencies for the sticker-spacer architecture could suppress the coalescence of multiple metastable droplets into a single large droplet.

Exhaustion of specific interaction valencies is a root cause of arrested LLPS even in the absence of active processes

First, we studied self-assembly by multivalent polymers whose adhesive domains interact via stable, ‘non-transient’ interactions to understand the early events in the growth process. We observed that, except for extremely high concentrations where polymers form a system spanning network (Figure 2A), the most feasible scenario at smaller concentrations is that of co-existing clusters with a concentration-dependent distribution of cluster sizes. While the non-monotonic dependence of density-transitions on concentration is consistent with equilibrium predictions (Harmon et al., 2017), the prevalence of multiple clusters with a distribution of cluster sizes at lower Cmono is a key deviation. Crucially, in the parameter regime where we observe dense clusters (Figure 2), the available interaction valencies get consumed within smaller-sized assemblies (Figure 3), making them inert for further cluster growth. Our results suggest that two critical time-scales decide whether a cluster continues to grow further – a) concentration-dependent time-scale of two chains (or clusters) encountering each other and forming the initial functional interaction, and, b) the exhaustion of valencies within a small cluster, a time-scale dependent on intrinsic features of the polymer (Figure 4). Crucially, these two time-scales are sensitive to subtle modifications in the self-assembling polymer chain (Figure 5).Therefore, while the underlying equilibrium landscape drives the process of self-assembly, the metastable multi-droplet system could be a dynamic outcome that is characteristic of polymers with multiple stickers. Our findings are consistent with the predictions from the ’sticky reptation’ theory proposed by Semenov and Rubinstein that suggests that the dynamics of associative polymers with multiple stickers exhibits a slow-down for high degrees of association (Rubinstein and Semenov, 2001). At higher degrees of association for multivalent polymers, with fewer free sites, the effective lifetime of each interaction becomes higher, making the recruitment of a new chain progressively less likely (Rubinstein and Semenov, 2001; Leibler et al., 1991). The role of physical crosslinks has also been reported to play a key role in tuning the viscoelasticity of protein assemblies (Roberts et al., 2018).

Modifications to the linker region can also result in altered densities of the self-assembled state, with a 10–100 fold enrichment in molecular concentrations within the clusters (Figure 5 and Figure 5—figure supplement 3), consistent with the experimentally reported degrees of enrichment within condensates (Nott et al., 2015; Mitrea and Kriwacki, 2016; Li et al., 2012) and previous theoretical studies showing the role of the linker in modulating phase separation (Harmon et al., 2017). Using LD simulations involving reversible and irreversible functional interactions, we establish a potential physical mechanism determining microphase separation in membrane-less organelles, with the finite nature of the specific interactions driving the peculiar phenomenon. Therefore, the highly cross-linked nature of the phase-separated state (Harmon et al., 2017; Rubinstein and Semenov, 2001) would result in the early droplets not being able to transition into a single large droplet even for reversible, transient interactions.

Bridging the gap between the early and biologically relevant timescales

To reach the time scales relevant to biology, we employed a coarse-grained kinetic model where each polymer (from the LD simulations) is represented as a diffusing reaction centre on a 2D lattice, which can interact either non-specifically (mimicking inter-linker interactions) or specifically with neighboring centers. The difference between the two types of interactions is that the number of specific interactions that each centre can make is limited by its valency. In an extension of the LD model, specific interactions are stable yet reversible and can form and break with rates dictated by detailed balance. In contrast to the 2D-lattice aggregation models employed by Dine et al., 2018, our model incorporates interaction valencies explicitly. Our kMC simulations also incorporate diffusion of entire clusters to allow coalescence, unlike the Dine model where clusters can grow only by monomer diffusion events and Ostwald ripening – an extremely slow process. Consistent with our LD simulations, our kinetic Monte-Carlo simulations of the phenomenological model reveal that for time-scales relevant to biology, fully percolated, system-spanning clusters (LclusNtot—> Ntot) are observed only for high bulk density (ϕlattice ) of monomers (Figure 7). Further, the phenomenon of exhausted adhesive valencies is more prominent for species with lower valency (fewer adhesive domains in the prototypical polymer), as evident from much smaller sizes of the largest cluster after an hour-long (real time) simulation run. The sensitivity of the cluster sizes to these phase parameters (valency, bulk density, interaction strengths) suggests that sizes and morphology of membrane-less organelles are amenable to tight control. In the context of the size-dependent function of membrane-less organelles, this tunability becomes critical. This is consistent with previously proposed mechanisms of organelle growth control by limiting the availability of components critical for droplet growth (Goehring and Hyman, 2012).

The kinetically trapped nature of the multi-droplet phase was further confirmed by the observation that cluster sizes (for the same simulation time) decrease upon slowing down the diffusion-limited time-scale by varying solvent viscosity (LD) or phenomenological diffusion coefficient (kMC) (Figure 9). The lattice-diffusion model, despite its minimalistic nature, captures the experimentally established relationship between molecular valency and critical concentration (Li et al., 2012) for phase separation (Figure 7B). An interplay between the valency of the generalized polymer and the strength of interactions can also alter the exchange times of molecules with the bulk medium dramatically. A slight shift in either these valencies or interaction strengths could result in a change in exchange rates by orders of magnitude. Further, for regions of the parameter space (ϵsp and λ) that favor large clusters (SclusNtot —> Ntot) separation, we observe a dramatic slow-down in molecular exchange times. In other words, for parameters that result in a fast-exchanging condensed state, the metastable microphase is the most favored outcome (Figure 8B and Figure 8—figure supplement 1). Such a discontinuity makes these systems extremely sensitive to mutations (Wang et al., 2018) that might cause a shift in dynamics and eventual loss of function of these droplets. Also, such shifts could also make these systems extremely responsive to non-equilibrium processes such as RNA-processing, side-chain modifications (acetylation, methylation) that are often attributed to modulating condensate dynamics (Wang et al., 2018; Hofweber and Dormann, 2019; Wang et al., 2014). The differential exchange times in response to variation of interaction parameters in our model lend further support to the scaffold-client model (Xing et al., 2018). The scaffolds that are slower exchange species with higher valencies acting could, therefore, recruit faster exchange clients with lower valencies. The valencies and strength of interactions could have thus evolved to achieve exchange times that ensure the functionality of spatial segregation via liquid-liquid phase separation.

A multi-domain sticker-spacer architecture allows for separation of two functions with the folded domains (conserved) performing functional role while the spacer regions being modified over time to tune the propensity to phase separate and also the material nature of the condensate. Overall, our multi-scale study shows that the block co-polymer-like organization of these multivalent proteins, with finite specific interactions driving phase separation, could manifest itself in the metastable multi-droplet state (Figure 10). A switch in the driving force for self-assembly, from the specific to non-specific interactions via sticky linkers (Wang et al., 2018), could alter not only the kinetics of assembly but also have implications in disorders associated with aberrant phase separation (Ramaswami et al., 2013; Patel et al., 2015). Our study paves the way towards the rational design of phase-separating polymers which is of vital importance for addressing their role in disease and function.

Graphical Summary.

Muli-valent proteins can exist in two Flory-like equilibrium states, the fully mixed solvated state and a single large macro-phase separated state (when free monomer concentration (Cmono) exceeds the critical concentration (Csat)). Above the critical concentration, the proteins first phase separate into multiple droplet phases (black arrow) which are long-living, metastable structures. However, the growth of these individual droplets to a single large macro-droplet (blue arrow) through coalescence and Ostwald ripening is unlikely at biologically relevant timescales. The multiple metastable micro-droplets can also age into highly cross-linked structures (green arrow) that saturate interactions within themselves and thus cannot grow further. Extremely stable functional interactions or sticky inter-linker interactions could also result in these arrested microphases.

Testable predictions

The emphasis of the current study is on the metastable nature of the long-living micro-droplet-like state of multivalent biopolymers. A signature of metastability would be the dependence of the observable quantities such as droplet sizes and size distribution on the dynamics of self-assembly process – a phenomenon not observed at equilibrium. Here, we propose a set of potential experiments that could be performed to test our key prediction of the dynamic origin of the multi-droplet system.

  • Our LD and kMC simulations suggest that the observed droplet size distributions (for a fixed biologically relevant time-scale) would depend on the viscosity of the solvent. At higher viscosities, the slowing down of droplet coalescence would therefore result in a shift in the droplet size distribution towards smaller sizes. A systematic experimental study of phase separation (at a time-scale of hours) in a range of solvent viscosities would provide an essential test of our findings.

  • Our phenomenological kinetic model predicts that the dynamic, droplet phase exists in a very narrow parameter regime (concentration, valency, interaction strength) and that a subtle change in parameters could result in either a fully mixed-phase or clusters that exhibit a dramatic slowdown in exchange times. A more systematic experimental study of droplet dynamics as a function of interaction valencies would be required to establish the veracity of these predictions. Engineered multivalent proteins such as the SH3-PRM system (Li et al., 2012) could serve as useful model systems for such a study.

  • Another key observation of this model is the ability of the linker region to modulate the nature of phase separation. Previously, synthetic multivalent protein model systems have been used to study LLPS in vitro (Li et al., 2012). Given that our dynamic model predicts small, condensed aggregates for strong inter-linker interactions, systematic studies of droplet size distributions in a range of inter-linker interaction strengths would establish how the observable quantities are sensitive to intra- and inter-cluster dynamics.

Materials and methods

Langevin Dynamics simulations

Force field

Request a detailed protocol

The polymer chains in the box are modelled using the following interactions. Adjacent beads on the polymer chain are connected via harmonic springs through the potential function,

(1) Estretching=ksi=1M1(|riri+1|r0)2,

ri and ri+1 refer to the adjacent ith and (i+1)th bead positions, respectively. Here, r0 is the equilibrium bond length and ks represents the spring constant. This interaction ensures the connectivity between the beads of a polymer chain.

This spring constant, ks , in our LD simulations is set to 5 kT/Å2, in order to ensure rigid bonds between neighboring beads of the same chain. The equilibrium length r0 , for bonds connecting adjacent linker beads, is set to 4.5 Å while the same for specific interactions between functional domains is set at 20 Å.

To model semi-flexibility of the polymer chain, any two neighboring bonds within the linker regions of the polymer chains interact via a simple cosine bending potential

(2) Ebending=κi=1M2(1cosθi),

where θi describes the angle between ith and (i+1)th bond while κ is the energetic cost for bending (bending stiffness of linker). We vary the bending energy parameter (κ) to model rigid (κ>2 kcal/mol) or flexible linkers (κ=2 kcal/mol).

The non-bonded isotropic interactions between linker beads and linker-functional interactions were modelled using the Lennard-Jones (LJ) potential,

(3) Enb=4ϵnsi<j[(σ|rirj|)12(σ|rirj|)6],

for all |ri-rj|<rc < rc, where rc refers to the cutoff distance beyond which the non-bonded potentials are neglected. The LJ potentials were truncated at a distance cutoff of 2.5σ. σ was set to 4.2 Å (roughly that of amino acid) for linker beads and 20 Å for functional domains.

The strength of the attractive component of this potential, ϵns, was varied to achieve varying degrees of inter-linker interactions in our simulations. In our simulations, the linker regions participate in inter-linker interactions only. The strength of this inter-linker ϵns in our simulations was varied in the range of typical strengths of short-range interactions in biomolecules. We vary the ϵns for pair-wise inter-linker interaction (per pair) in the range of 0.1 kcal/mol (≈0.2 kT) to 0.5 kcal/mol (≈1 kT) (Sheu et al., 2003). For comparison, strong non-covalent interactions such as the H-bonds are known to be of the order of 0.5 kcal/mol to 1.5 kcal/mol in solvated proteins. Similar values of isotropic, short-range interactions have been used in conventional coarse-grained protein force-fields such as the MARTINI (Marrink et al., 2007; Monticelli et al., 2008).

Modeling specific interactions

The specificity of the functional interactions is modelled by imposing a valency of 1 between different complementary specific bead types (red and yellow beads in Figure 1) via a bonding vector attached to each bead. Valency of 1 per site means that each specific interaction site can only be involved in one such interaction. The total valency of an individual polymer chain (λ) is the number of adhesive sites that belong to a single chain. Bond formation (modelled using stochastically forming harmonic springs) occurs with a probability (Pform) if two complementary beads are within a defined interaction cutoff distance (rc). The spring constant for functional interactions is set at 2 kT/Å2. In simulations with irreversible functional interactions, a bond, once formed, remains stable for the duration of the simulations. In order to implement reversible specific interactions, we introduce a bond breakage probability (Pbreak) when the inter-domain distance stretches beyond a breakage cutoff – rbreak. In our reversible specific interaction simulations, a bond breaks with a probability (Pbreak ) of 1 when the two domains move 2.2 Å from the equilibrium distance of the spring (used to model the specific interaction). The strength of the spring at the breakage point is ≈4 kT, a value that is within the range of interaction strengths previously probed in lattice simulations of sticker-spacer polymers (Harmon et al., 2017).

LD simulations details

Request a detailed protocol

The dynamics of these coarse-grained polymers was simulated using the LAMMPS molecular dynamics package (Plimpton, 1995). In these simulations, the simulator solves for the Newtons’s equations of motion in presence of a viscous drag (modeling effect of the solvent) and a Langevin thermostat (modeling random collisions by solvent particles) (Schneider and Stoll, 1978). The simulations were performed under the NVT ensemble, at a a temperature of 310 K. The mass of linker beads was set to 110 Da while the mass of the idealized functional domains (red and yellow beads in Figure 1) was set to 7000 Da that is approximately equal to the mass of the SH3 domain. The size of the linker beads was set at 4.2 Å (of the same order as amino acids) while that of the functional domains was set at 20 Å (≈size of a folded SH3 domain Musacchio et al., 1992). The viscous drag was implemented via the damping coefficient, γ=m/6πηa. Here, m is the mass of an individual bead, ‘η’ is the dynamic viscosity of water and ‘a’ is the size of the bead. An integration time step of 15 fs was used in our simulations, and the viscosity of the surrounding medium was set at the viscosity of water (10-3 Pa.s). Similar values of these parameters have been previously employed for coarse-grained Langevin dynamics simulations of proteins (Bellesia and Shea, 2009).

Kinetic Monte carlo simulations

To assess biologically relevant time-scales, we develop a phenomenological kinetic model wherein the individual multivalent polymer chains are modelled as diffusing particles with fixed valencies (Figure 6). In order to study the temporal evolution of a system of diffusing particles on a lattice, we employ the (Gillespie, 1977) or the kinetic Monte Carlo simulation, a technique used for modelling a system of chemical reactions. This algorithm has previously been used to model biological processes as diverse as gene regulation (Parmar et al., 2014) and cytoskeletal filament growth (Ranjith et al., 2009).

Reactions (processes), and interactions modelled in the simulation

Request a detailed protocol

We begin the simulation with a set of Ntot particles randomly placed on a 2D square lattice. The density of particles on the lattice (ϕlattice ) is a measure of the bulk initial concentration of monomers in these simulations. Each particle in our lattice Monte carlo simulations is a coarse-grained representation of the bead-spring polymer chains in the LD simulations, with an effective valency (λ) that is a simulation variable. Assuming that there are unoccupied neighboring sites, the monomers can diffuse on the lattice in any of the four directions (in 2D) with a rate kdiff. Particles occupying adjacent sites on the lattice experience a weak, non-specific interaction (of strength ϵns) that is isotropic in nature. This attractive isotropic force is analogous to the inter-linker interactions in the LD simulations. A functional bond (specific interaction) can stochastically form between any pair of neighboring particles, provided both particles possess unsatisfied valencies. The rate of bond formation between a pair of particles with unsatisfied valencies is kbond. On the other hand, an existing functional bond can break with a rate kbreak that is equal to kbond.exp(−ϵsp ) in magnitude. ϵsp is the strength of each functional interaction. Additionally, the entire cluster that any given monomer is a part of can diffuse in either of the four directions with a scaled diffusion rate that is inversely proportional to the size of the cluster (kdiff /Sclus ). A particle that is part of a cluster can diffuse away from the cluster with a rate kdiff.exp(-ϵsp + ϵns ). Here, ϵsp + ϵns is the magnitude of net interactions that any particle is involved in. These rates are schematically described in Figure 6.

Details of the Gillespie simulation algorithm

Request a detailed protocol

At every time-step of the simulation, these rates – diffusion of monomer, diffusion of the cluster that a monomer is a part of, rate of bond formation and rate of bond breakage – are initialized for every monomer based on the current configuration of the system. Note that, for any given configuration, all or only a subset of these events could be possible. We then advance the state of the system by executing one reaction at a time. The probability of each event is computed, and it is equal to ri/rij . Here i and j refer to the identity of the monomer and the event type (diffusion, cluster diffusion, bond formation) , respectively. The denominator in this term is the sum of the rates of all the possible events for a given configuration (for all monomers), referred to as rtotal. At any given time-step, given the set of probabilities for each event, we choose the event to be executed by drawing a uniformly distributed random number z1 (between 0 and 1). The configuration of the system is then updated based on the event that gets chosen at any time-step. The Gillespie Exact Stochastic Simulation technique (Gillespie, 1977) is a variable time-step simulation approach wherein the simulation time is advanced using the following expression, δt=(1/rtotal)ln(z2), where z2 is a uniform random number. Here, rtotal is the sum of the rates of all possible events at a given simulation time-step. This δt is based on the assumption that the waiting times between any two events are exponentially distributed. The above algorithm has to be iterated several times such that each reaction has been fired multiple times, suggesting the system has reached steady-state behaviour. A detailed account of the algorithm can be found in the seminal work by Gillespie, 1977.

Computing cluster sizes in kMC and LD simulations

Cluster size calculations form the basis of the current study. In order to compute cluster sizes in LD simulations, we used the gmx-clustsize module within the GROMACS molecular dynamics package (Van Der Spoel et al., 2005). When the closest distance between any two functional domains from different polymer chains is within a 22.5 Å radius (interaction radius for functional domains in the simulation), we consider them to be part of the same cluster. We use the snapshots from the final 1 μs of the trajectory (and from five independent trajectories) to perform cluster size computations. In order to compute cluster sizes and distributions in the kMC simulations, we employed the Hoshen-Kopelman algorithm (Hoshen and Kopelman, 1976), which is routinely used for studying percolation problems.

Rationale behind phase parameters in KMC simulations

Request a detailed protocol

Two key phase parameters in our kMC simulations are valency per interacting particle and bulk density of particles on the lattice. We performed simulations for valencies ranging from 3 to 6, typical of multivalent proteins forming membrane-less organelles. Unlike non-specific interactions which can be only one per neighbor, thereby a maximum of 4 per particle on a 2D-lattice, there could be more than one specific interactions between a pair of neighboring particles, as long as both participating members have unsatisfied valencies. This multiplicity of specific interactions between nearest neighbors in this model reflects the fact that each polymer contains λ>1 globular domain that can engage in specific interactions with globular domains from neighboring protein. The bulk density of particles is defined as, ϕlattice = Ntot/L2, where Ntot and L are the number of particles and lattice size, respectively). In our kMC simulations, we varied ϕlattice within the range of 0.01 to 0.1, a range that is consistent with the analogous parameter for LD simulations, the occupied volume fractions within the LD simulation box. For instance, a free monomer concentration (Cmono) of 10 μM corresponds to a volume fraction of 0.008, which increases to 0.17 for 200 μM. Further, the phase diagrams for KMC simulations were primarily computed for a weak non-specific ϵns of 0.35 kT, chosen in order to focus on specific interaction-driven cluster growth. In the kMC simulations, ϵns refers to the net non-specific interaction for a pair of monomers as opposed to a pair of interacting linker beads in the LD simulations. It must, however, be noted that the interaction strength for a pair of interacting chains is not merely additive (with the length of the polymer chain). A potential explanation for this could be found in earlier studies wherein it has been argued that smaller segments (≈ length of 7–10 amino acids) within a long, solvated polymer chain (like IDPs) could behave like independent units referred to as blobs (Pappu et al., 2008). Since the degree of coarse-graining employed in kMC is of the order of one particle per polymer chain, it lacks the microscopic degrees of freedom of the bead-spring polymer chain (in LD simulations). Hence we do not employ a higher non-specific interaction strength in our phase plot computations. As a proof of principle, we show the mean pairwise interaction energies for dimers from the LD simulations, for a weak inter-linker interaction strength of 0.1 kcal/mol (Figure 7—figure supplement 4A) and a corresponding phase-diagram for non-specifically driven interactions (Figure 7—figure supplement 4B). The mean cluster sizes and molecular exchange times in the kMC simulations were computed over 100 independent kMC trajectories.


  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55

Decision letter

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

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

There is growing interest in and a cognitive dissonance between the equilibrium expectation of a single large droplet vs. multiple micro-phase separated droplets that one observes from phase separation that gives rise to biomolecular condensates. Several explanations have been put forth to explain extant observations. Your work provides an intriguing and testable hypothesis based on the dynamics of phase separation of multivalent molecules that conform to a stickers-and-spacers architecture. The concept of two distinct timescales namely, the timescale of diffusion that brings molecules into contact and the reconfiguration times before which the available valence is exhausted provides another explanation for the observation of micro-phase separated droplets as opposed to a single large droplet.

Decision letter after peer review:

Thank you for submitting your article "Membraneless organelles are dynamic, metastable long-living droplets formed by sticker-spacer proteins" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Rohit V Pappu as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by José Faraldo-Gómez as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Alexander Grosberg (Reviewer #3).

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

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, when editors judge that a submitted work as a whole belongs in eLife but that some conclusions require a modest amount of additional new data, as they do with your paper, we are asking that the manuscript be revised to either limit claims to those supported by data in hand, or to explicitly state that the relevant conclusions require additional supporting data.

There is growing interest in the topic of "liquid-liquid phase separation" better described as phase separation aided bond percolation as a mechanism by which multivalent macromolecules form biomolecular condensates. The sequence-intrinsic features that give rise to biomolecular condensates can be measured via in vitro reconstitutions and instantiated in numerical simulations. The general framework whereby multivalent macromolecules are mapped onto a stickers-and-spacers architecture is proving to be quite useful for describing the phenomenology and thermodynamic driving forces for phase transitions specifically., phase separation aided bond percolation. Despite these advances, the actual dynamics of phase separation aided bond percolation have been ignored to this point. Are there interesting insights that emerge from a detailed understanding of the dynamics of phase separation aided bond percolation? The answer, according to Ranganathan and Shakhnovich appears to be a resounding yes! Through a combination of multiscale computations, the authors show that multivalent proteins can become dynamically arrested in so-called "micro-phases". As a result, instead of the equilibrium picture of a single, gigantic condensate coexisting with a dispersed dilute phase, one observes the formation of numerous miniature condensates coexisting with the dilute phase. These micro-phases come about because of the rapid manner in which the available valence is exhausted when compared to the slower timescales associated with clusters of molecules diffusing, colliding and coalescing. This gives rise to a dynamical phase diagram that deviates from the equilibrium picture, although the equilibrium free energy landscape still underwrites the driving forces as well as the energetic and / or entropic barriers associated with the dynamics of phase separation aided bond percolation. The results presented by the authors are likely to garner interest given the apparent cognitive dissonance between the equilibrium picture that can be recapitulated in vitro and the observation of multiple smaller facsimiles of condensates coexisting with one other in vivo. This problem has engendered active and animated discussions in the field, and the authors' results as well as their perspectives offer new food for thought.

The manuscript, in its original form, has been reviewed separately by three reviewers. The reviews have been discussed and the general consensus is that this is an interesting and timely contribution that builds on recent successes of the stickers-and-spacers formalism for describing equilibrium phase behavior of multivalent macromolecules. The reviewers noted the importance of the observations that emerge from the simulations. Reviewer 1 is concerned about some bold assertions being made that appear to invalidate the equilibrium picture. It is worth emphasizing that the dynamics are underwritten by the intrinsic free energy landscape and in this work, which does not include sources or sinks, that should be it. That there are interesting dynamics that cause deviations from the equilibrium picture is something that has been reported in a completely different context (please see: https://www.nature.com/articles/s41563-018-0182-6). The consensus is that the authors should articulate the distinction between the equilibrium phase diagram and the dynamical phase diagram, highlight the differences that emerge from the two, and strive to minimize sweeping statements regarding the (in)validity of one over the other, since one (the equilibrium picture) begets the other (the dynamical picture). Reviewer 1 highlighted the crucial need for analyzing density transitions (ignored in the current work) as well as networking transitions (as has been done by Harmon et al. and Choi et al.). This is going to be crucial to understand if the arrested phases are in a distinct regime as far as density is concerned, even though they do appear to be above the percolation threshold as far as connectivity is concerned. The reviewers were unanimous in their dual concerns regarding the number of jargon terms that were introduced throughout the narrative, the interchangeable usage of these jargon terms making things look imprecise, and the assertive use of terms like gel, macroscopic gel, solid-like phases with no quantification of how one can make such assertions absent rheological data (see comments by reviewer 3), characterization of density transitions (see comments by reviewer 1), or experimental data (see comments by reviewer 2). There were concerns about the inaccessibility of the methods descriptions that were raised by reviewer 2 with regard to the description of the Langevin dynamics simulations as well as the Kinetic Monte Carlo simulations. An underlying concern, articulated primarily by reviewer 2, pertains to the artificial imposition of irreversible crosslinks in the LD simulations. The implied concern is that by imposing irreversibility, the observations rather become a tautology. Indeed, while this will not affect the percolation threshold calculated in the mean-field limit, it almost certainly will impact the true estimate of the percolation threshold. A justification for the approach and preferably a comparison to a more realistic scenario in the LD simulations is warranted. Reviewer 2 expressed concerns regarding the hard to follow description of the KMC simulations. A series of concerns were raised and these should be addressed via revisions. Finally, reviewer 3 raised concerns regarding the description of reptation and its usage as a concept in analyzing the simulation results. This is crucial because it sets up the separation of time scales and reviewer 3 is of the opinion that a) the formalism as captured in the mathematics is erroneous and b) an alternate picture might be more relevant to describe the data. Reviewer 1 also highlights the fact that concept of arrested phases is not new. Lee and colleagues have written a series of papers highlighting the possible role of active processes as determinants of micro-phases in cells. The work of Boeynaems et al. highlights the formation of arrested phases due to stable RNA structures. These efforts should be discussed.

As for the title, the assertive declaration is misleading. This work does not actually pertain to a specific biomolecular condensate or membraneless organelle. The title should be modified to better capture the fact that a computational model of multivalent macromolecules modeled using the stickers-and-spacers framework leads to the prediction of long-lived micro-phases. This is the essence of the MS and the title should reflect this point precisely.


All in all, while there is strong support for seeing this work published in eLife, there are numerous concerns that could be addressed in a revised version with new analysis sans new data, although the inclusion of reversible crosslinks in the LD simulations would be very helpful, even essential. There is concern that the findings will be opaque to the average reader of eLife. Therefore, it would be helpful if the authors were to lean on a few colleagues who are primarily molecular or cellular biologists and gain their perspective regarding the overall accessibility of the revised MS and ensure that the message will be clearly understood.

Essential revisions:

The major revisions requested are a) please moderate the rather extreme pronouncements about the irrelevance of the equilibrium picture; b) please analyze density as well as percolation transitions; c) please replace the usage of gelation with bond percolation; d) please streamline the narrative to use a small number of necessary jargon terms and define these up front; e) please redo the LD simulations with reversible crosslinks to obtain a comparative assessment of the pictures that emerge with irreversible vs. reversible crosslinks; f) please reconsider the assertions being made about solid-phases since the results cannot be used to make these assertions; g) please note that Ostwald ripening has been observed in vitro and there are numerous examples, especially from the Brangwynne and Hyman labs showing that a large single droplet will form and coexist with a dispersed dilute phase; h) please include citations to the work of Chiu Fan Lee, Boeynaems et al. and Roberts et al. https://www.nature.com/articles/s41563-018-0182-6; i) please revisit the description and usage of the reptation model, correct the mathematical description to ensure that what emerges has units of time, and ensure that usage of this model is indeed appropriate / accurate; j) please provide a more accessible treatment and an illustrative example of the KMC formalism (details that enable reproducibility are important); k) please provide a summary of the main conclusions, please minimize some of the overstatements and casting aside of previous work, cite the seminal work of Semenov & Rubinstein, and touch base with the work of Huan-Xiang Zhou, especially the recent publication in PNAS.

Reviewer #1:

In this work, the Ranganathan and Shakhnovich present their analysis of dynamical aspects of phase transitions for associative polymers described by a stickers-and-spacers architecture. This framework has been recently adapted to describe the phase behavior of linear and branched multivalent proteins. The central findings in the current work are that multivalent molecules form dynamically arrested micro-phases and that this comes about because of network terminating interactions arising at short time scales that in turn inhibit the growth of the network to give rise to macro-phase separation wherein a single condensate (droplet) defined by a percolated network within the condensate coexists with a dispersed dilute phase. Instead of this equilibrium scenario, the disparate timescales lead to multiple droplets that act as metastable sinks to inhibit the growth of macroscopic phases. In the picture that emerges, the dynamics of phase transitions of stickers-and-spacers systems are governed by a hierarchy of timescales and the separation between distinct timescales gives rise to stable mesophases i.e., metastable phases. This aspect of the authors' finding is a) interesting, b) timely and warrants reporting, and c) reasonably well established in the synthetic polymer literature, although under-appreciated (perhaps even unappreciated) in the biomolecular condensate literature. This work certainly deserves to be published in eLife. There are, however, a few key issues that would benefit from careful scrutiny.

1) The authors suggest that the equilibrium state is rarely accessed in vivo or in vitro. This flies in the face of data that have been presented in vitro. Observations of Ostwald ripening and numerous reports of a single large droplet forming in test tubes that coexist with droplet-free dispersed phases have been reported in the literature. There are technical difficulties associated with observing / recording growth dynamics in vitro, especially when the assays are based on fluorescence microscopy, but careful analysis from several labs show that the equilibrium scenario is indeed realized in vitro. One recent example of this is the following paper: https://science.sciencemag.org/content/367/6478/694, but I emphasize that there are many more. What is however true (and not mentioned in the current MS) is that there have been several discussions in the literature about apparent emulsification being observed in vivo. The work of Lee, Jülicher and coworkers introduced the idea of active emulsions – see https://iopscience.iop.org/article/10.1088/1361-6633/ab052b for a review, a direct use of this theory for modeling stress granules http://hdl.handle.net/10044/1/57849, and the original work that introduced the idea of active processes giving rise to essentially uniform distributions of droplets https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.078102. Simply stated, there is an active (no pun intended) debate about the mechanisms by which microphases arise. The authors have introduced a dynamical view of this could arise and they suggest that this could be general, which may be the case. However, it is worth emphasizing for accuracy and rigor that the presence / absence of arrested phases will be governed by a combination of factors including what the authors articulate here as well as the quench depth into the two-phase regime (which requires prior knowledge of the equilibrium phase diagram), and the rigidity effects due to prior structure formation – see for example https://www.pnas.org/content/early/2019/03/28/1821038116. The upshot is that it is incorrect to assert that extant data support the contention that the equilibrium scenario is rare nor is it accurate to suggest that the equilibrium phase behavior has zero bearing on "reality". In fact, the equilibrium phase diagram provides the crucial touchstone (the intrinsic free energy landscape) for the dynamical phase diagram. In light of these facts, it would be useful to provide a more nuanced representation of the importance of the equilibrium framework and spell out what the dynamical phase diagram adds to the field (which is a lot). One needs both descriptions, not either / or.

2) In the manuscript, the term phase separation is used often, but the work does not ever analyze an order parameter for density transitions, focusing instead on the clustering alone. This is important because the narrative ends up conflating networking / clustering (gelation i.e., bond percolation), which is what is being quantified, with density transitions i.e., phase separation, which is never actually quantified. The works of Harmon et al., 2017, and Choi et al., 2019, go to great lengths to distinguish between phase separation and bond percolation. In the current context, this would seem to be especially relevant. Are the microphases clusters with densities that are in accord with the dilute phase or do the densities also change? If the latter, then how do the densities (concentrations) within the dilute phases compare to the thermodynamic saturation concentration and how do the densities within the microphases compare with the equilibrium densities within dense phases? These questions can only be answered by deploying two sets of order parameters, one based on networking (which the authors use) and one based on density transitions (ignored in this work).

3) Please note that the phase diagrams presented here are dynamical phase diagrams. They are not equilibrium phase diagrams and to suggest otherwise is misleading. This should be fixed.

4) Finally, on a semantic note: the manuscript is quite long, takes a rather narrow view of the literature, and importantly, introduces too many over-loaded terms (terms used differently in different contexts) and too many distinct terms as well. It would help immensely to streamline the verbiage to define a set of terms and stick with these throughout. Otherwise, I am concerned that the average reader of eLife will be quite confused and not dial into the central message of this work.

Reviewer #2:

In this work, the authors use a computational approach to model phase separation of polymers that are comprised of stickers and spacers (similar to that employed in Harmon et al., 2017). By employing computational Langevin dynamics and Monte Carlo simulations, the authors monitor phase separation on biologically-relevant timescales. In so doing, the authors are able to explore the different tunable parameters of phase separation, including non-specific interactions among spacer regions (linkers), specific interactions among sticker regions, linker rigidity, and polymer concentration. By adjusting these parameters, their simulations revealed several interesting points. First, as expected, large clusters (consisting of many polymers) form only at high monomer concentrations. Second, phase separation occurred only when there was high valency with functional interactions among stickers, and not in the absence of such interactions. Third, simulations showed "arrested droplets", i.e. no system-spanning macro-phase separated state was observed at intermediate monomer concentrations. The authors provide evidence that this is due to exhausted free valencies. This observation is tunable based on the interaction strength of the linkers and stickers. The LD simulations also allow for linker rigidity modeling that informed on timescales that impact the timing of kinetically-arrested droplets. Importantly, the authors identify competing timescales for how clusters grow, specifically (1) the initial diffusion encounter, and (2) the reorganization time for the clusters. Next, the authors employed kinetic Monte Carlo simulations to more accurately model interactions that take place. Above, when polymers encountered each other, interactions were irreversible, i.e. once two polymers interacted, they did not come apart again for the rest of the simulation. This is a poor assumption, therefore, the second half of the work employs reversible binding. However, this dramatically changes how the simulation models polymer dynamics, as the polymer is now represented coarsely as a single bead with a variable λ parameter that designates occupied valency.

Overall, the manuscript provides some key insights into the kinetics of droplet assembly, as well as into the material properties of droplets as a function of time. This study further illustrates the usefulness of the stickers and spacers polymer model to obtain the rules for how polymers self-assemble into droplets at relevant timescales. The manuscript needs to be strengthened by clearer language, a better presentation of the results from varying many different parameters that could be tested experimentally, careful summary of the results, and detailed methods so that the general audience can appreciate how the conclusions were reached.

1) There is a need to pictorally represent the “metastable” simulation states when the valencies are filled. This is necessary to illustrate the difference between the macro-phase separated state of a system-spanning single droplet, and what is meant by micro-phase separation of many droplets. Care needs to be taken with the term phase separation propensity as this manuscript quantifies the parameters necessary for phase separation; why not use saturation concentration thresholds (csat), and discuss increasing/decreasing csat? These terms need to be clearly defined for the readership.

2) There is a lot of useful information presented throughout the manuscript of how various parameters tune the kinetics of droplet assembly in the simulation. It is necessary to sum these key points in a conclusion section that is missing from the manuscript (although this reviewer appreciates the testable hypotheses section for useful experimental follow-up experiments). That being said, it would have been very useful to have experimental data to back up some of the simulation results, perhaps via a collaboration with experimental groups on the SH3-PRM system.

3) There are many parts of the text that could be improved with clearer definitions. For example, subsection “Phenomenological kinetic simulations predict microphase separation at biologically relevant timescales” is hard to read, and needs to be broken up so that the results are described succinctly. A second example: I imagine there will be confusion regarding what is considered “gel” in Figure 7 for readers that are familiar with the term in polymer physics, vs. biologists that may confound the definition of “gel” with solid-like material properties.

4) There are interesting results in Figure 5 regarding the relationship between Rg and the concentration threshold for phase separation. These results are particularly interesting for evaluating the contribution of linker regions in multidomain proteins. Are these findings applicable for all phase separating systems? I started thinking about this in the context of proteins that phase separate with increasing temperature (LCST) vs. proteins that phase separate with decreasing temperature (UCST).

5) There is difficulty in reading the Materials and methods section for the kinetic Monte Carlo simulations, and more discussion in the text should be devoted to explaining how the method works, particularly for a general audience. Perhaps an example timestep would be useful here. It is unclear to me how the simulation time is advanced. What is meant by rtotal? Are these the sums of individual rates across all monomers and event types for a given time?

6) Along these lines, please carefully reread the LD methods section. Multiple variables are missing from the text (following Equation 4). This made it impossible for me to understand how LD simulations were executed, and I would like to see this section again.

7) Implications from kinetic MC simulations are that the interplay between interaction strength of specific interactions and valency tune the liquid-like properties of phase separated droplets. Increasing specific interaction strength promotes more solid-like droplets. This is an important point, given recent studies that discuss effects of lysine/arginine residues on droplet material properties and csat. Can the authors comment on this point in their manuscript, even though the sticker-spacer model does not specifically address different amino acid types? This would add value to the findings.

8) Figure 8 needs example simulation results to help the reader understand u1 and u2 differences.

9) Scripts detailing the analysis for the simulations (e.g. largest cluster analysis) are missing and need to be included in the manuscript for clear reproducibility.

Reviewer #3:

This is an interesting paper dealing with a currently very hot subject of liquid-liquid phase separation inside the cell with the participation of appropriate proteins. The heart of the work is a Langevin dynamics simulation accompanied and illustrated by a simple semi-analytical model. The most important novel results deal with the prediction of microphase separation (unlike classical Flory theory predictions) due to the temporal exhaustion of available valences. Importantly, authors claim that this type of phase segregation happens on a biologically relevant time scale.

In my opinion, the main results seem novel, sound, and physically transparent. The model is a reasonable compromise between tractability and approach to reality. I would definitely recommend the work for publication. This general conclusion notwithstanding, I have two issues with the work.

1) One is the use of some loose terms, such as "gel-like" phase; what exactly is it? Even more, there is also "dynamically solid macroscopic gel-like state" as opposed to "liquid-like". What do all these words mean? As far as I understand, authors did not measure shear modulus of their aggregates, so their characterization in terms of solid versus liquid requires clarification. Furthermore, if gel is a percolating covalent network, then it is strictly speaking a solid, at least on the scale above its mesh size. I would welcome a more concise usage of all these words.

2) My second remark is about reptation. Estimate of the time called reptation is very important for the general conclusions of this work. This estimate comes apparently from formula (1) which is very strange. First of all, units do not seem to work in this formula: L = 145 is the number of beads, it is dimensionless, while Dbead is the diffusion coefficient; the result cannot be time. While this may seem like a simple typo (e.g., forgotten bead size squared), combined with other issues this one may not be that simple. Second, the real reptation time, unlike formula (1), is not controlled by Dbead but by an L times smaller quantity. Third, perhaps the strangest, this strange formula is supported by the reference Feldman, 1989, which is actually a review of the famous textbook by Doi and Edwards. The book would be of course the most welcome source about reptation, but not a review of this book. In fact, my suspicion is that the process in question is not reptation, but some sort of a diffusion along the chain. In any case, authors should really clarify what exactly they mean.

To summarize, I would feel comfortable recommending this work for publication in eLife as soon as authors clarify the two issues mentioned above.


Author response

Essential revisions:

The major revisions requested are:

a) Please moderate the rather extreme pronouncements about the irrelevance of the equilibrium picture.

We appreciate the concerns raised by the reviewers with regards to the lack of acknowledgment of the importance of equilibrium theories. In the revised manuscript, at several places in the Introduction as well as the Discussion, we highlight how the process of phase-separation by multi-valent polymers is driven by the underlying thermodynamic landscape. In the revised manuscript, we emphasize that the current study focuses on potential role of dynamics in hindering the progression of a multi-droplet system to the equilibrium two-state system. We also provide examples from literature supporting the observation of droplet coalescence in support of the equilibrium predictions. The current study sheds light on the importance of dynamics of cluster growth, complementing the existing equilibrium understanding of the phenomenon.

b) Please analyze density as well as percolation transitions.

We thank the reviewer for this suggestion. In the revised manuscript, we introduce a second order parameter to analyze the intracluster density of polymers (normalized by the bulk density). This quantity is analogous to the order parameter “ρ” used to analyze density transitions by Harmon et al1. In Figure 1A of the revised manuscript, we plot the intracluster densities and cluster sizes as a function of total monomer concentration (Cmono). As observed in case of equilibrium lattice simulations by Harmon et al1, the density transition shows a non-monotonic behavior as a function of concentration, with large concentrations resulting in system-spanning networks with low densities. In a narrow range of concentrations, we do observe dense clusters. Further, in Figure 5, we show how intracluster density depends on the properties of the linker. In Figure 5B and Figure 5—figure supplement 2, we show how inter-linker interactions can tune the density of polymers within the clusters.

c) Please replace the usage of gelation with bond percolation.

We thank the reviewers for this suggestion. In the revised manuscript, we refer to large clusters with system-spanning networks as a “macrophase” (characterized by low density within the macro cluster), and to a system of coexisting clusters (Sclus << Ntot) as a metastable “micro-phase”. We refrain from using the term gelation in the revised manuscript.

We define these terms in the first subsection of the Results titled “Terminology and Notations used in this study”.

d) Please streamline the narrative to use a small number of necessary jargon terms and define these up front.

We now introduce the conventions and terminology used in the article in the first subsection of the Results titled “Terminology and Notations used in this study”. We also introduce a table listing the physical interpretation and the definition of all variables and order parameters used in the current study.

e) Please redo the LD simulations with reversible crosslinks to obtain a comparative assessment of the pictures that emerge with irreversible vs. reversible crosslinks.

To test whether the early time-scale behavior changes in the presence of reversible interactions, we introduced breakable functional bonds in our model. A comparative assessment of the single largest cluster sizes as well as the distributions for reversible and irreversible crosslinks is now presented in Figure 2B and C. As with the irreversible interaction simulations, even in the presence of breakable interactions, Lclus << Ntot (Figure 2B) indicating the existence of long-living metastable microphase droplets except for large Cmono when we observe a system-spanning network. Critically, we find a coexistence of intermediate cluster sizes with small and large clusters suggesting an increased diversity in cluster sizes at the early stages of droplet assembly even upon introduction of breakable interactions (Figure 2C). However, a more detailed study of impact of bond formation and breakage dynamics on cluster growth is the subject of future work.

The implementation of reversible functional interactions is explained in the “Modeling Specific Interactions” subsection under Materials and methods.

f) Please reconsider the assertions being made about solid-phases since the results cannot be used to make these assertions.

The usage of “solid-like” and “liquid-like” phases in the original manuscript was in reference to the rate at which monomers get exchanged between the condensate and the free medium in our kMC simulations (referred to as exchange times in the manuscript). We do not perform any other detailed assessment of the material properties of these assemblies. We therefore do not use these terms in the revised manuscript. Instead, we call these refer to these phases as slow- or fast-exchange phases in the revised manuscript , based on their mean exchange times in kMC simulations.

g) Please note that Ostwald ripening has been observed in vitro and there are numerous examples, especially from the Brangwynne and Hyman labs showing that a large single droplet will form and coexist with a dispersed dilute phase.

Thanks for pointing this out. In the Introduction section of the revised manuscript, we now include references supporting the in vitro observation of droplet growth via Ostwald ripening and coalescence.

h) Please include citations to the work of Chiu Fan Lee, Boeynaems et al. and Roberts et al. https://www.nature.com/articles/s41563-018-0182-6

In the Introduction and Discussion sections, we now broaden the discussion to include the potential role of active mechanisms as well as the rigidity of the scaffold in promoting/stabilizing a multi-droplet system. We cite the suggested references in the relevant context within the revised text.

i) Please revisit the description and usage of the reptation model, correct the mathematical description to ensure that what emerges has units of time, and ensure that usage of this model is indeed appropriate / accurate.

We agree with the reviewer’s comments about the validity of the reptation based mathematical description of the dimer reorganization timescale. Indeed, a classical reptation model predicts scaling of N3 while our formula predicts Rouse-type scaling of N2. Extra power of N in reptation model comes from slowing down diffusion in entangled dense polymer solution (melt) which does not play a role in our estimate of time scale of valency saturation in a polymer dimer in solution. However, we feel that this interesting question warrants a separate study. We removed this equation and accompanying discussion from revised manuscript since it is not central to our message here. We plan to return in the near future to careful study of polymer physics of valency exhaustion in dilute solution and report findings in a separate manuscript.

We also cite the sticky-reptation model described by Semenov and Rubinstein in the Discussion section of the revised version. While relevant, it explores somewhat different model of a dense solution while our focus here is on dynamic processes leading to formation of metastable microdroplets from initially dilute solution.

j) Please provide a more accessible treatment and an illustrative example of the KMC formalism (details that enable reproducibility are important).

We now provide additional details of the kMC simulation, and the details of cluster size computations in the Materials and methods section of the revised manuscript. The reactions modeled in the kMC simulations are schematically described in Figure 6. The simulation code for the lattice-kMC simulations will also be deposited in a software repository and will be accessible to the readers.

k) Please provide a summary of the main conclusions, please minimize some of the overstatements and casting aside of previous work, cite the seminal work of Semenov & Rubinstein, and touch base with the work of Huan-Xiang Zhou, especially the recent publication in PNAS.

The summary of the main conclusions is presented in the form of a succinct illustration (Figure 10) as well as the “Testable Predictions” section. We believe that a separate additional conclusions section will lengthen the paper even further and duplicate already existing pictorial conclusion. The work by Semenov and Rubinstein is now discussed in the context of exhausted valencies within clusters in the Discussion section. We also discuss the patchy particle model by Huan-Xiang Zhou’s group in the Introduction section while referring to existing computational studies related to LLPS. We thank the reviewers for bringing these important works to our attention.

Overall, we modified three figures in the main text (Figure 2, Figure 4 and Figure 5) and included a new subsection introducing the key terminology and order parameters used in this study. We include a new subsection in the Materials and methods section, discussing the implementation of breakable interactions. Under Introduction and Discussion in the revised manuscript, we also provide a more elaborate context emphasizing the importance of the existing thermodynamic understanding. Also, to provide a broader view of existing literature, we incorporate references to the possible role of active mechanisms in shaping droplet size distributions and compare predictions from both mechanisms.


Article and author information

Author details

  1. Srivastav Ranganathan

    Department of Chemistry and Chemical Biology, Harvard University, Cambridge, United States
    Conceptualization, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
  2. Eugene I Shakhnovich

    Department of Chemistry and Chemical Biology, Harvard University, Cambridge, United States
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4769-2265


National Institutes of Health (RO1 GM068670)

  • Eugene I Shakhnovich

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


This work was supported by NIH RO1 GM068670. We are grateful to Alexander Grosberg, Mark Miller and Ranjith Padinhateeri for useful discussions.

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


  1. Rohit V Pappu, Washington University in St Louis, United States
  2. Alexander Grosberg, New York University, United States

Publication history

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


© 2020, Ranganathan and Shakhnovich

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.


  • 529
    Page views
  • 128
  • 0

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)

  1. Further reading

Further reading

    1. Cell Biology
    2. Physics of Living Systems
    Anjali Jaiman, Mukund Thattai
    Research Article Updated
    1. Immunology and Inflammation
    2. Physics of Living Systems
    Marco Molari et al.
    Research Article Updated