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

Figure 1

Download asset Open asset

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.

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.

When $L_{clus}\to N_{tot}$—> N_tot, we refer to the state as a system-spanning ‘macro-phase’. A system of multiple, co-existing clusters with with $S_{clus}\ll N_{tot}$<< N_tot, 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 (${\varphi }_{clus}$ /$\varphi$) 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.

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 $C_{mono}$, 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 $C_{mono}$(Figure 2—figure supplement 1A, purple curve) . Therefore, in this regime of ${ϵ}_{ns}$ and $C_{mono}$, 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 ($C_{mono}$). The polymer-chains, for the data plotted in Figure 2, have flexible linkers ($\kappa$=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 ($C_{mono}$<70 μM), $S_{clus}\ll N_{tot}$<< N_tot. 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 – ${\varphi }_{clus}$ /$\varphi$. 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 $C_{mono}$, 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 ${\displaystyle {\varphi }_{clus}/\varphi \to 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

Download asset Open asset

Figure 2—source data 1

Compressed zip file containing the source data for cluster sizes and cluster size distributions (along with raw unprocessed data) plotted in Figure 2.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig2-data1-v2.zip

Download elife-56159-fig2-data1-v2.zip

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, $L_{clus}\ll N_{tot}$ << N_tot (Figure 2B) except for large $C_{mono}$ 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.

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 $S_{clus}\ll N_{tot}$ << N_tot, 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 $C_{mono}$ (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, $P_{form}$ , 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

Download asset Open asset

Figure 3—source data 1

Compressed zip file containing the source data for temporal evolution of different species (dimer,trimer,monomer etc) for different concentrations plotted in Figure 3.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig3-data1-v2.zip

Download elife-56159-fig3-data1-v2.zip

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

Figure 4

Download asset Open asset

Figure 4—source data 1

Compressed zip file containing the source data for Figure 4.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig4-data1-v2.zip

Download elife-56159-fig4-data1-v2.zip

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 ($\kappa$<2 kcal/mol in Figure 4B), the linkers behave like flexible polymers. In this limit ($\kappa$<2 kcal/mol), the mean time for the valency to exhaust within the dimer (referred to as the re-organization time, $T_{exhaust-valency}$ ) is independent of linker stiffness. For ($\kappa$>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 $T_{exhaust-valency}$ with $\kappa$. 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 $T_{exhaust-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.

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 $C_{mono}$ 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 (<$R_g^{chain}$ >) 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 <$R_g^{chain}$ >, 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 $C_{mono}$ 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 $C_{mono}$ 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

Download asset Open asset

Figure 5—source data 1

Compressed zip file containing the source data for largest cluster size, density of the cluster and radius of gyration of constituent chains, for different values of ${ϵ}_{ns}$ plotted in Figure 5.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig5-data1-v2.zip

Download elife-56159-fig5-data1-v2.zip

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 ($T_{exhaust-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 $T_{exhaust-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 $T_{exhaust-valency}$. Our dimerization simulations show that a second mechanism of slowing down the $T_{exhaust-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.

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 $N_{tot}$ such multivalent particles (at varying densities, ${\varphi }_{lattice}$ ) that diffuse freely, at a rate $k_{diff}$, and the number of particle collisions per unit time is proportional to $k_{diff}$*${\varphi }_{lattice}$ . In our kMC study, we vary ${\varphi }_{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 $k_{bond}$. The valency per particle ($\lambda$) 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 $k_{break}$ = $k_{bond}$ *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 ($k_{diff}$ /$S_{clus}$). 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, $k_{diff}$ and $k_{bond}$. The second time-scale, $T_{exhaust-valency}$, is a time-scale that depends on the $k_{bond}$ and $k_{break}$ 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.

Figure 6

Download asset Open asset

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 (${\varphi }_{lattice}$), (b) rate of bond formation ($k_{bond}$), (c) valency per interacting particle ($\lambda$), and (d) the strength of specific interactions (${ϵ}_{sp}$). With the assumption that, for diffusion-limited self-assembly, the rate of free diffusion $k_{diff}$ is the fundamental time-scale limiting cluster growth, we first explore the relationship between the rate of specific bond formation ($k_{bond}$), and $k_{diff}$ (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 ${\displaystyle k_{bond}/k_{diff}\to 0}$, 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, $L_{clus}\ll$<<1) for low-intermediate ${\varphi }_{lattice}$ even at biologically relevant time-scales of hours in the kMC simulations. In this regime of low-intermediate ${\varphi }_{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 ${\varphi }_{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

Download asset Open asset

Figure 7—source data 1

Compressed zip file containing the source data for kinetic phase diagrams in Figure 7.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig7-data1-v2.zip

Download elife-56159-fig7-data1-v2.zip

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 ${\displaystyle k_{bond}/k_{diff}\to 0}$). For comparison, in Figure 7—figure supplement 4, we present the mean cluster sizes for assemblies that are stabilized by non-specific interactions only ($\lambda$=0). The ${ϵ}_{ns}$-${\varphi }_{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 ($k_{bond}\to k_{diff}$) establishing the validity of the instantaneous bonding assumption in the LD simulations. It is the ratio of $k_{bond}$/$k_{diff}$ , and not the absolute magnitudes, that is a vital parameter for these simulations. Hence, in all the kMC simulations we set the value of $k_{diff}$ to 1 $s^{-1}$ and vary the ratio $k_{bond}$/$k_{diff}$ 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 $\lambda$ and ${\varphi }_{lattice}$ as the phase parameters. For smaller $\lambda$ and low ${\varphi }_{lattice}$ , $L_{clus}$ (and $S_{clus}$ ) $\ll N_{tot}$N_tot (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 $\lambda$ 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 ($\lambda$=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}$ -$k_{bond}$ 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.

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 $\lambda$=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 $\lambda$=3 suggesting that an interplay between $\lambda$ and ${ϵ}_{sp}$ could tune the droplets for desired exchange properties. We further varied these two parameters ($\lambda$ 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 $\lambda$ 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 $\lambda$ 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

Download asset Open asset

Figure 8—source data 1

Compressed zip file containing the source data for kinetic phase diagrams in Figure 8.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig8-data1-v2.zip

Download elife-56159-fig8-data1-v2.zip

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 ($\eta$=${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 $\lambda$ 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

Download asset Open asset

Figure 9—source data 1

Compressed zip file containing the source data for mean largest cluster sizes (and size distributions) from LD and MC simulations studying the effect of solvent viscosity.

https://cdn.elifesciences.org/articles/56159/elife-56159-fig9-data1-v2.zip

Download elife-56159-fig9-data1-v2.zip

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.

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 $C_{mono}$ 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.

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 ($L_{clus}\to N_{tot}$—> N_tot) are observed only for high bulk density (${\varphi }_{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 $\lambda$) that favor large clusters ($S_{clus}\to N_{tot}$ —> N_tot) 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.

Figure 10

Download asset Open asset

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.

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 ($\lambda$) is the number of adhesive sites that belong to a single chain. Bond formation (modelled using stochastically forming harmonic springs) occurs with a probability ($P_{form}$) if two complementary beads are within a defined interaction cutoff distance ($r_c$). The spring constant for functional interactions is set at 2 kT/Ų. 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 (P_break) when the inter-domain distance stretches beyond a breakage cutoff – $r_{break}$. In our reversible specific interaction simulations, a bond breaks with a probability ($P_{break}$ ) 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).

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

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.

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.

© 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.