Introduction

To fulfill their biological functions, cells must organize their contents into different compartments. One way of achieving spatiotemporal organization is via the formation of membrane-bound organelles. Another, that does not rely on membranes, is the self-assembly of proteins and nucleic acids into biomolecular condensates through the process of liquid–liquid phase separation (LLPS)1,2,3. In LLPS, a protein solution demixes into a protein-rich liquid phase—the condensate—in coexistence with a protein-poor liquid matrix. Condensates are ubiquitous within cells, with some of the most famous examples including stress granules4, P granules5, nuclear speckles6, and the nucleolus7.

The interest in understanding the molecular grammar8, 9 and biophysical determinants of cellular LLPS has significantly increased in recent years due, in part, to the realization that condensates take part in a wide-range of cellular functions—including genome silencing10, 11, signaling12, buffering cellular noise13, and formation of super-enhancers14, among many others15, 16. Furthermore, neurodegenerative diseases such as Alzheimer’s, Parkinson’s or ALS have been associated with aberrant or misregulated phase transitions17,18,19.

The chief drivers of biomolecular LLPS are multivalent molecules (i.e., those possessing multiple interaction sites) known as ‘scaffolds’20. Scaffolds are defined as the molecules—e.g., a single type of protein, a combination of various cognate proteins, or mixtures of proteins and nucleic acids21—that are able to sustain LLPS on their own. Above a critical concentration, scaffolds establish a sufficiently high number of transient attractive scaffold–scaffold interactions per unit of volume (the ‘liquid-network connectivity’22) to compensate for the entropic loss upon demixing. In addition, scaffolds can establish interactions with biomolecules that are unable to phase separate on their own, known as ‘clients’, recruiting them to the condensates via scaffold–client interactions. The recruitment of clients can alter significantly the stability and structural properties of condensates22,23,24. Therefore, together, scaffolds and clients shape the biophysical properties of the biomolecular condensates that play fundamental roles within the cell25,26,27,28.

A striking observation is the presence of condensates inside cells that do not grow beyond a certain size7, 29,30,31. Basic thermodynamics suggest that over time, LLPS should result in the formation of a single large condensate rather than multiple coexisting small droplets32. The latter case is disfavored because, cummulatively, it yields a high surface area to volume ratio (i.e., high interfacial free energy penalty), while the former ensures that the interfacial free energy penalty is minimized. Despite of the thermodynamic preference for single condensates over multidroplet systems, both types of architectures are present within cells. Indeed, single large condensates have been observed26, 33, but also diverse coexisting size-restricted droplets have been reported in different in vivo systems, such as in the amphibian oocyte nucleolus7, in the Saccharomyces cerevisiae cytoplasmic processing bodies34 and elsewhere35,36,37. Additionally, size-conserved multidroplet architectures have been found in non-coalescing ribonucleoprotein condensates26 and in multiphase complex coacervates in vitro38. However, the underlying molecular mechanisms and biophysical driving forces behind the formation of multiple coexisting droplets, also known as emulsification39, require further investigation.

Several explanations for why or how emulsions can be thermodynamically stable at biological relevant time-scales are currently under debate40. One possible mechanism is that the presence of active ATP-dependent processes might conveniently regulate the conditions where droplets grow and coalesce41,42,43,44. Other studies suggest that proteins with various highly distinct interacting domains may form micelle-like condensates31, 45, 46. In multicomponent mixtures, another possibility could be that specific binding proteins act as powerful surfactants and, thus, reduce the droplet surface tension penalty leading to multicondensate coexistence35, 36, 38, 47. Moreover, a recent alternative explanation suggests that the interplay between protein diffusion and saturation of protein binding sites can also induce size-conservation in condensates29, 30. It is plausible that all these different mechanisms contribute to the size conservation of condensates inside cells under different conditions.

In this work, we use a minimal protein model, which recapitulates the experimentally-observed relationship between protein valency and critical parameters23, 48, 49, to investigate the regulation of droplet size in binary mixtures of multivalent proteins (scaffolds and clients). We show that liquid–liquid phase separation of scaffolds and clients mixtures, where clients act as surfactants, can give rise to single droplets or multiple size-conserved droplets (Fig. 1a). Further, we reveal that the transition between the two scenarios can be regulated by the condensate scaffold/surfactant client ratio. Our simulations suggest how general molecular features such as protein connectivity, binding affinity, and droplet composition can critically modulate and stabilize the formation of size-conserved condensates29, 30, and might have also implications to understanding the phase behavior of multilayered condensates23, 38, 50.

Results

A minimal protein model for scaffold and client mixtures

Coarse-grained potentials have emerged as powerful tools for describing the phase behavior of biomolecules, such as proteins and nucleic acids, and delineating the underlying physicochemical features that drive LLPS51, 52. Various levels of molecular resolution can be achieved with coarse-grained models; encompassing mean field descriptions53, 54, lattice-based approaches55, 56, minimal models23, 48, 49, 56,57,58,60, and sequence-dependent representations46, 60,61,62,63,65. Here, we employ our minimal protein model48, which has been previously applied to unveil the role of protein multivalency in multicomponent condensates49 and multilayered condensate organization66, as well as to investigate the role of RNA in RNA-binding protein nucleation and stability67. In this model, proteins are described by a pseudo hard-sphere potential68 that accounts for their excluded volume, and by short-range potentials for modeling the different protein binding sites, and thereby mimicking protein multivalency48 (Fig. 1b). For computational efficiency, an implicit solvent is used; therefore, the condensed phase corresponds to a liquid phase, and the diluted phase to a vapor. In what follows, the unit of distance is \(\sigma\), the molecular diameter of the proteins (both scaffolds and clients), and the unit of energy \(k_{\text B}T\) (for further details on the model parameters and the employed reduced units see the “Methods” section).

Following the framework of Banani et al.20, scaffolds are defined as proteins that can establish both homotypic interactions and heterotypic interactions with clients, while clients (or hereafter called surfactants) as proteins that are limited to bind only to scaffolds (i.e., they do not bind to other clients; except where otherwise stated). Within this scheme, phase separation is driven by scaffolds (high-valency proteins), whereas surfactants (proteins with lower valency) are recruited to condensates at the expense of depleting LLPS-stabilizing scaffold–scaffold interactions. We model scaffolds as 4-binding sites particles and surfactants as 3-binding sites particles. (Fig. 1b). This choice fulfills two important requirements: (1) it ensures that scaffolds have a higher valency than the clients, and (2) it allows us to easily establish a common simulation temperature for the system at which each type of patchy particle distinctly behaves as either a scaffold or a surfactant client. The two-phase coexisting densities as a function of temperature for our model of scaffold proteins (blue), a 50:50 binary mixture of scaffolds and surfactants (black), and a system composed of just surfactants (red) are depicted in Fig. 1c (Top panel). Note that while in the 50:50 mixture, surfactant proteins do not interact homotypically, in the pure surfactant system they do; the latter allows us to compare the effect of protein valency on the critical parameters of pure self-interacting protein systems. These simulations show that the addition of surfactant proteins that are strong competitors for the scaffold binding sites significantly hinders the ability of scaffolds to phase separate (i.e., clients lower the critical temperature)20, 49. Moreover, the presence of surfactant clients drastically reduces the surface tension of the condensates (black curve) as shown in Fig. 1c (Bottom panel). In the following section, we elucidate the implications of the client-induced surface tension reduction on the behavior of phase-separated condensates.

Figure 1
figure 1

(a) Schematic representation of single-droplet formation versus size-conserved protein condensation. Scaffold proteins are depicted as blue spheres, while surfactant clients are shown as red spheres. (b) Minimal coarse-grained model for protein LLPS: Blue and red spheres represent the excluded volume of scaffold and surfactant (client) proteins respectively, while gray patches represent the binding sites of the proteins. Two different proteins are modeled: scaffold proteins, with 4 promiscuous binding sites in a tetrahedral arrangement, and surfactant proteins, with 3 binding sites in a planar equidistant arrangement that can only bind to scaffold binding sites (except where otherwise stated). Details on the model parameters are provided in the “Methods” section. (c) Top: Phase diagram in the temperature–density plane for a scaffold protein system (blue), a 50:50 scaffold–surfactant mixture (black) and for a hypothetical surfactant system in which client proteins can self-interact (red). The same self-interacting potential employed for scaffold proteins (blue) is also applied for surfactant proteins that can hypothetically self-interact (red). This system serves to further illustrate the effect of multivalency in LLPS. Filled circles indicate the estimated coexisting densities from Direct Coexistence simulations68,69,71, and empty ones depict the critical points calculated via the law of rectilinear diameters and critical exponents72 (see Supporting Information for details on these calculations). Bottom: Surface tension (\(\gamma\)) dependence on temperature for the scaffold system (blue), the 50:50 scaffold–surfactant mixture (black), and the hypothetical system in which surfactant proteins can self interact (red). Filled squares represent direct estimations of \(\gamma\), continuous curves depict fits to our data of the following form: \(\gamma \propto (T^*-T^*_{\text {c}})^{1.26}\), and empty squares show the critical temperature of each system evaluated through the law of rectilinear diameters and critical exponents72. The vertical dashed line indicates the temperature at which the remainder of our study is performed. Note that temperature (in reduced units, \(T^*\)) is renormalized by the critical temperature \(T^*_{\text {c}}\) (also in reduced units) of the scaffold protein system (\(T^*_{\text {c}}=0.12\)).

Full size image

Size-conserved condensates can be regulated by surfactant proteins

Despite their seemingly subordinate role in condensate formation, client molecules (e.g. low valency proteins or self-avoiding nucleic acids) can significantly impact the organization of biomolecular condensates20, 38, 73, 74. Remarkable examples include the multilayered organization of stress granules75, 76, the nucleoli26, nuclear speckles77, mixtures of RNA-binding proteins with RNA78, 79, and in vitro complex coacervates38, 73, 74. While scaffolds mainly contribute to maximizing the condensate liquid network connectivity20, 49, clients can severely reduce the droplet interfacial free energy cost of creating an interface23. The proliferation marker protein Ki-67, a component of the mitotic chromosome periphery, is an important example of a natural surfactant in intracellular compartmentalization35. In this context, our simulations show that low-valency clients that do not establish homotypic interactions behave as surfactants within condensates, as they can greatly decrease the condensate surface tension [Fig. 1c (Bottom)].

To further investigate the role of surfactant clients, we perform two sets of Direct Coexistence (DC) simulations68,69,71: (1) for the pure scaffold system, and (2) for the 50:50 binary mixture of scaffold and surfactants shown in Fig. 1c. At a constant temperature of \(T^*/T^*_{\text {c}}=0.75\) and a global (system) density of \(\rho ^*=0.136\) (both in reduced units; see “Methods” section for further details on the employed reduced magnitudes), we create three different simulation box geometries with the dimensions summarized in Table 1. Using this approach, we can effectively modulate the surface/volume ratio of the condensed phase (hereafter called droplet): where the droplet surface is \(S = 2*L_y*L_z\) and its volume is \(V = L_y*L_z*L_{x,slab}\); with \(L_{x,slab}\) representing the width of the condensate in the x direction, and \(L_x\), \(L_y\) and \(L_z\) being the different sides of the simulation box. Figure 2a–c, summarizes the phase behavior of the pure scaffold system (Top panels) and the 50:50 binary mixture (Middle panels) along the three designed simulation box geometries. The pure scaffold condensate exhibits distinct surface/volume ratios depending on the box geometry (see Table 1), while maintaining the same droplet density on all three cases. In other words, the scaffold condensate would be able to continuously grow as a single droplet at expense of the diluted phase until reaching equilibrium. In contrast, the scaffold–surfactant mixture yields various coexisting equilibrium droplets with a roughly constant surface/volume ratio of 0.21(2) \(\sigma ^{-1}\) in all systems (Table 1). This result strongly suggests how size-conserved multidroplet phase behavior can be simply induced by the presence of low valency clients, which actively contribute to lowering the surface tension of the droplets [Fig. 1c (Bottom)]—thereby serving as natural surfactants35, 36.

Table 1 Simulation box dimensions and condensate surface/volume ratios (S/V) of the three box geometries represented in Fig. 2 for the pure scaffold system and the 50:50 binary mixture of scaffold and surfactant proteins. Geometries (a), (b) and (c) account for the (a), (b) and (c) panels shown in Fig. 2, respectively.
Full size table

We next analyze the composition of the different coexisting scaffold–surfactant condensates along the distinct box geometries [Fig. 2 (Bottom panels)]. In all cases, including the pure scaffold condensate, the properties of the droplets are remarkably similar; i.e., the density of all droplets, as well as their composition and surface tension are roughly constant (Table 2). Notably, we find that the surfactant density profile becomes higher than that of scaffolds at the droplet interface, showing how the partition coefficient of surfactants is greater than that of scaffolds in the outer region. Previous works suggest that accumulation of surfactants at the interface is preferable as it minimizes the condensate surface tension35, 38, 66. We also verify that the presence of multiple coexisting droplets is the thermodynamically stable state, rather than just metastable, by simulating over sufficiently long timescales to allow for multiple droplet fusion events and variations in droplet composition. Importantly, these tests reveal that even when in contact the droplets coexist without coalescing or altering their equilibrium composition. The multidroplet behavior of size-conserved condensates (in our case of \(\sim\) 0.21 \(\sigma ^{-1}\) surface/volume ratio) is a consequence of the thermodynamic conditions of our system (i.e., mixture composition, temperature, and density). Note that droplet curvature effects such as Laplace internal pressure80 or surface tension dependence on droplet curvature81 have not been considered in our simulations, since we do not expect them to play an important role at biologically relevant droplet size scales (\({\mathscr {O}} \mu\)m)40. Those effects are only expected to be dominant in the nanometer scale (i.e., up to droplet radii of tens of nanometers)81,82,84.

The presence of surfactant clients within the condensate substantially lowers the liquid network connectivity49, 85, 86 and, therefore, reduces the enthalpic gain sustaining LLPS. Consequently, the system minimizes its free energy by optimizing the number of surfactants that are incorporated into the condensed phase; i.e., by creating higher surface/volume ratios, where surfactants are preferentially located towards the interface rather than in the core. Such free energy optimization yields multiple coexisting condensates of a certain size, rather than a single-condensate system. Moreover, the emergence of multiple coexisting droplets stabilised by surfactant proteins instead of a single condensate can also contribute to increase the entropy of the system, and thus, lowering its free energy. Nonetheless, this behavior is only thermodynamically favorable when the condensate surface tension is very low, as in the case of the 50:50 scaffold–surfactant condensates with a surface tension of \(\gamma =0.23\,k_{\text B}T\,\sigma ^{-2}\)—almost 7 times lower than that of the pure scaffold condensate (\(\gamma =1.58\,k_{\text B}T\,\sigma ^{-2}\)) at the same conditions (\(T^*/T_c^*=0.75\)).

Figure 2
figure 2

Direct Coexistence simulations for a scaffold protein system (Top panels) and a 50:50 binary mixture of scaffold and surfactant proteins (Middle panels) with different simulation box geometries (see Table 1 for the different simulation box sides employed in a, b and c geometries). All simulations were performed at \(T^*/T^*_{\text {c}}=0.75\) and a global density of \(\rho ^*=0.136\). (a)–(c) (Top and Middle panels): Scaffold proteins are colored in blue and surfactant clients in red; the protein binding sites of both protein types are colored in gray. (a)–(c) (Bottom panels): Density profiles along the long side (x) of each simulation box Lx for the 50:50 scaffold–surfactant mixture. The (scaffold+surfactant) total density profile is depicted in black, the surfactant client density profile is shown in red, and that of scaffolds in blue. Densities of the different scaffold–surfactant coexisting droplets, as well as the scaffold molar fraction in the condensed and diluted phase, for each geometry are given in Table 2. The density of the pure scaffold condensates in all geometries is 0.54(1).

Full size image
Table 2 Properties of systems presented in Fig. 2 containing 50:50 scaffold–surfactant compositions in different simulation box geometries. In all cases, 6 independent simulations (with different initial velocity distributions), each starting from a pre-equilibrated configuration, were performed. For the geometries with more than one droplet, the values are averaged over the different droplets, although the variance between distinct droplets is significantly small as shown in Fig. 2 (Bottom panels).
Full size table

Surfactant concentration critically modulates droplet size

As discussed in the previous section, our minimal protein model shows that for a given composition of scaffold and surfactant clients, independent of the imposed box geometry, droplets can only grow to a certain size. This size-restricted growth is, in turn, determined by the optimal surface/volume ratio that minimizes the free energy of the coexisting liquid phases. Therefore, larger system sizes lead to higher number of coexisting size-conserved droplets. On the other hand, when the condensate is only composed of scaffold proteins, as the system size increases, the size of the condensate simply grows instead of yielding new multiple size-conserved droplets (Fig. 2). These results illustrate how a simple model for scaffold and surfactant proteins, merely controlled by protein valency and binding affinity, can recapitulate mesoscale features of in vivo and in vitro condensates that exhibit size-conserved growth7, 26, 34,35,36,37,38. Since this phase behavior only arises when the concentration of surfactants is not negligible, we now investigate how condensates can switch between both scenarios, and how their surface/volume ratio is modulated by their relative scaffold–surfactant composition.

Figure 3
figure 3

Surface tension (\(\gamma\)) dependence on the scaffold–surfactant composition (in \(\%\)) at \(T^*/T^*_{\text {c}}=0.75\) and system (global) density of \(\rho ^*=0.136\). Vertical dashed lines indicate the maximum surfactant concentration which allows LLPS for a given number-droplet regime in our system. Note that the maximum droplet size varies continuously with surfactant concentration even within the same number-droplet regime. The maximum condensate size in terms of surface/volume ratio along the different number-droplet regimes are: \(\sim\) 0.11 \(\sigma ^{-1}\) for the one-droplet regime, \(\sim\) 0.19 \(\sigma ^{-1}\) for the two-droplet regime, and \(\sim\) 0.27 \(\sigma ^{-1}\) for the three-droplet regime. At surfactant client compositions exceeding 65\(\%\) there is no LLPS. For each computed composition, the global density of the system and the simulation box cross-section (\(L_z L_y\)) are kept constant. Snapshots of the DC simulations along different surfactant composition droplet regimes are included. Please note that in the multidroplet regime, the different coexisting droplets (at the same surfactant concentration) exhibit similar compositions and densities.

Full size image

By gradually increasing the surfactant client concentration of the scaffold–surfactant mixture (at a constant temperature and system density), the condensate size progressively decreases to accommodate the equilibrium droplet surface/volume ratio to the simulation box geometry; i.e., the condensate splits into two and, subsequently, into three coexisting liquid droplets (see Fig. 3). We note that the composition of the different coexisting droplets at a given surfactant concentration is remarkably similar, highlighting that all droplets are in equilibrium. In parallel, we evaluate the surface tension of the different coexisting droplets as a function of surfactant concentration. We find that \(\gamma\) monotonically decreases (but not linearly) as the surfactant concentration increases (Fig. 3). This result is not surprising (see Fig. 1c) given that one of the key molecular driving forces behind size-conserved multidroplet formation is the reduction of \(\gamma\) by surfactants coating the droplet surface. Above a certain surfactant concentration—exceeding 65 % for our system and at the given conditions—LLPS is inhibited. Beyond this limit, the condensate liquid network connectivity sustained by scaffold proteins can no longer compensate the mixing entropy of the system. We also analyze the surface/volume ratio of the droplets as a function of surfactant concentration. At infinitely low surfactant concentration, the condensate displays a ratio of \(\sim\) 0.09 \(\sigma ^{-1}\). However, such a ratio is fully determined by the total number of proteins in the system and the box geometry, since as shown in Fig. 2, scaffold condensates can reach any droplet size when surfactants are absent. At low surfactant concentrations (i.e., \(\%\) surfactant \(< 27.5\%\)), the maximum droplet size corresponds to a surface/volume ratio of \(\sim\) 0.11 \(\sigma ^{-1}\). Beyond that threshold concentration, the condensate shrinks, and to achieve the equilibrium surface/volume ratio, it splits into smaller coexisting droplets. The maximum equilibrium droplet size in the two-droplet regime is that corresponding to ratios of \(\sim\) 0.19 \(\sigma ^{-1}\). Finally, for surfactant compositions higher than \(38\%\), three coexisting condensates emerge. The maximum surface/volume ratio that droplets can achieve is \(\sim\) 0.27 \(\sigma ^{-1}\) at \(60\%\) client composition, which is only possible due to the extreme reduction in the surface tension (\(\gamma =0.15\,k_{\text B}T\,\sigma ^{-2}\))—more than one order of magnitude lower than that of the pure scaffold condensate (\(\gamma =1.58\,k_{\text B}T\,\sigma ^{-2}\)) at the same temperature and system density.

Figure 4
figure 4

(a) Droplet density profile (in reduced units) for scaffold (blue) and surfactant proteins (red) from the droplet center of mass towards the surrounding diluted phase for the 50:50 binary mixture within the simulation box geometry shown in Fig. 2a at \(T^*/T^*_{\text {c}}=0.75\). (b) Average number of engaged binding sites per protein (\(\varphi\)) as a function of distance from the droplet center of mass for scaffold (blue) and surfactant proteins (red). One binding site is considered to be engaged to another if the distance between them is less than 0.145 \(\sigma\) (i.e., the maximum bond length interaction between distinct protein binding sites; for further details on these calculations see Supporting Information).

Full size image

Previous studies have highlighted the challenges associated with measuring condensate surface tensions due to the small size of protein droplets3, 87. Nonetheless, there are available estimates of this magnitude for ribonucleoprotein condensates, and these measurements demonstrate that surfactant proteins can reduce \(\gamma\) by orders of magnitude26. With our minimal model, we qualitatively observe such behavior when surfactant proteins are recruited to the condensate, giving rise to emulsions of multiple coexisting droplets with very low surface tension. Surfactant proteins can lead to the formation of multidroplet emulsions by inducing multilayered condensate architectures26, 38. Diverse biomolecular organelles, either in the cell nucleus, such as the nucleoli26, or nuclear speckles77, as well as in the cytosol, such as stress granules,75, 76 exhibit this type of organization. Moreover, different in vitro complex coacervates38, 50, 73, 74, bioengineered ribonucleoprotein condensates in living cells88 and mixtures of RNA-binding proteins and RNA molecules78, 79 are also known to show multilayered assemblies. In Fig. 4a, we analyze the droplet architecture of a protein condensate with a 50:50 scaffold–surfactant composition [i.e., the system in Fig. 2a (Middle panel)]. We find that in the droplet core, the scaffold to surfactant ratio is much higher than along the interfacial region (Fig. 4a), where it drops to almost half that of the surfactant proteins. Nonetheless, the surfactant concentration within the droplet core is still remarkably high considering its destabilizing role in the condensate liquid network connectivity49. The observed non-homogeneous condensate organization stems from the higher valency of scaffold proteins, which allows them to establish higher molecular connectivity within the core condensate and, thus, induce higher enthalpic gain upon multilayered assembly.

Finally, to gain further insight on the droplet liquid network connectivity, we evaluate the average number of engaged binding sites per protein (\(\varphi\)) as a function of distance from the center of mass of the condensates (Fig. 4b). Scaffold proteins present a significantly higher amount of molecular connections per protein than surfactants (i.e., \(\varphi \sim 3.7\) and \(\varphi \sim 2.5\) for scaffolds and surfactants, respectively), at the droplet core. This observation highlights how surfactants negatively contribute to the stability of the condensate. However at the interface, such diminished connectivity of surfactant proteins (\(\varphi \sim 1\)) with respect to that of scaffolds (\(\varphi \sim 3\)) substantially reduces the condensate surface tension—by decreasing the enthalpic cost (\(\Delta h_{i}\)) of creating an interface. This energetically favorable protein arrangement, controlled in our model just by the variance in protein valency of the components, is expected to be contributed also by changes in relative binding affinities among proteins, and be modulated by post-translational modifications in vivo7, 34. Furthermore, such variations can be also relevant to understand the physical parameters controlling multilayered condensate organization and, ultimately, regulation of the formation of size-conserved multidroplet emulsions29,30,31.

Conclusions

In this work, we employ our minimal protein model48 to demonstrate how biomolecular multidroplet emulsification can be controlled by the subtle balance between liquid network connectivity and droplet surface tension, and how general molecular features such as protein valency and binding affinity can critically regulate this behavior. By using a binary mixture of scaffold and client proteins that act as surfactants (following the original definition proposed by Banani et al.20), we design a set of Direct Coexistence simulations in which we can conveniently modulate the simulation box geometry to assess the propensity of the condensates to accommodate different surface/volume ratios. The ability (or disability) of these mixtures to adopt different surface/volume ratios, imposed by the box geometry, can be regarded as an indirect measurement of the droplet propensity to grow beyond a certain size. We find that pure component scaffold condensates can easily adapt to distinct surface/volume ratios; in support of their ability to grow and fuse into a single droplet. However, 50:50 binary scaffold–surfactant mixtures stabilize instead several coexisting liquid condensates with roughly constant surface/volume ratios to accommodate the imposed system size geometry. Such behavior is a clear signature of size-conserved multidroplet emulsification, as found in the nucleolus7, ribonucleoprotein condensates26, 88, micelle-like condensates31, 45, and in vitro complex coacervates38.

We also elucidate the role of surfactant concentration in size-conserved droplet growth. By gradually decreasing the scaffold–surfactant ratio in our mixtures, we observe that the maximum droplet size is reduced, while simultaneously increasing the number of coexisting condensates. This trend continues until a sufficiently high surfactant concentration is reached, where LLPS is no longer possible. Moreover, as clients are added, the droplet surface tension dramatically decreases, facilitating the formation of multiple coexisting small liquid droplets at low interfacial energetic cost. Client proteins, besides decreasing the stability of the condensates49, can effectively behave as natural droplet surfactants35, 36. Due to their considerably lower molecular connectivity compared to that of scaffolds, surfactants preferentially migrate towards the droplet interface; thereby, minimizing the enthalpic cost of creating an interface23. Heterogeneous molecular organizations of condensates have been observed in stress granules75, 76, the nucleoli26 and nuclear speckles77. We find that such heterogeneity enable the maximization of the condensate liquid network connectivity through scaffold–scaffold protein interactions within the droplet core.

Rationalizing the underlying mechanisms employed by cells to precisely regulate the size of their diverse membraneless compartments and processing bodies7, 34, 36 represents a crucial step towards understanding intracellular spatiotemporal cell organization. Taken together, our coarse-grained simulations help to elucidate the relationship between single-droplet phase formation and size-conserved multidroplet architecture, and put forward general molecular features such as valency and binding affinity as chief drivers in these scenarios.

Methods

We model our coarse-grained multivalent proteins using the MD-Patchy potential proposed in Ref.48, which is composed by two different set of potentials: a Pseudo Hard-Sphere (PHS) potential89 to continuously describe the repulsive interaction and excluded volume between different protein replicas, and a continuous square-well (CSW)90 potential to describe the patch-patch interactions among different protein binding sites. The \(u_{PHS}\) potential is described by the following expression:

$$\begin{aligned} u_{PHS}= {\left\{ \begin{array}{ll} \lambda _r\left( \frac{\lambda _r}{\lambda _a}\right) ^{\lambda _a} \varepsilon _R \left[ \left( \frac{\sigma }{r}\right) ^{\lambda _r}- \left( \frac{\sigma }{r}\right) ^{\lambda _a} \right] + \varepsilon _R; &{}\quad \text{ if } r < \left( \frac{\lambda _r}{\lambda _a}\right) \sigma \\ 0; &{}\quad \text{ if } r \ge \left( \frac{\lambda _r}{\lambda _a}\right) \sigma \end{array}\right. } \end{aligned}$$
(1)

where \(\lambda _a=49\) and \(\lambda _r=50\) are the exponents of the attractive and repulsive terms respectively, \(\varepsilon _R\) accounts for the energy shift of the PHS interaction, \(\sigma\) is the molecular diameter (and our unit of length) and r is the center-to-center distance between different PHS particles. For the patch–patch interaction we use the following expression:

$$\begin{aligned} u_{CSW}=- \frac{1}{2} \varepsilon _{CSW} \left[ 1-tanh\left( \frac{r-r_w}{\alpha }\right) \right] \end{aligned}$$
(2)

where \(\varepsilon _{CSW}\) is the depth of the potential energy well, \(r_w\) the radius of the attractive well, and \(\alpha\) controls the steepness of the well. We choose \(\alpha =0.005\sigma\) and \(r_w=0.12\sigma\) so that each patch can only interact with another single patch.

The mass of each patch is a 5% of the central PHS particle mass, which is set to \(3.32 \times 10^{-26}\) kg, despite being this choice irrelevant for equilibrium simulations. This 5% ratio fixes the moment of inertia of the patchy particles (our minimal proteins). The molecular diameter of the proteins, both scaffold and clients, is \(\sigma =0.3405\) nm, and the value of \(\varepsilon _{R}/k_B\) is 119.81K. All our results are presented in reduced units: reduced temperature is defined as \(T^*=k_{\text B}T/\varepsilon _{CSW}\), reduced density as \(\rho ^*=(N/V)\sigma ^3\), reduced pressure as \(p^*=p\sigma ^3/(k_{\text B}T)\), and the reduced unit of time as \(\sqrt{\sigma ^2m/(k_{\text B}T)}\). In order to keep the PHS interaction as similar as possible to a pure HS interaction, we fix \(k_{\text B}T/\varepsilon _R\) at a value of 1.5 as suggested in Ref.89 (fixing \(T=179.71\)K). We then control the effective strength of the binding protein attraction by varying \(\varepsilon _{CSW}\) such that the reduced temperature, \(T^*=k_{\text B}T/\varepsilon _{CSW}\), is of the order of \({\mathscr {O}}(0.1)\).

Since both \(u_{PHS}\) and \(u_{CSW}\) potentials are continuous and differentiable, we perform all our simulations using the LAMMPS Molecular Dynamics package91. Periodic boundary conditions are used in the three directions of space. The timestep chosen for the Verlet integration of the equations of motion is \(\Delta t^*=3.7 \times 10^{-4}\). The cut-off radius of the interactions of both potentials is set to 1.17\(\sigma\). We use a Nosé–Hoover thermostat92, 93 for the NVT simulations with a relaxation time of 0.074 in reduced units. For NpT simulations, a Nosé–Hoover barostat is employed with the same relaxation time94.

The methodological details of the calculation of the phase diagram, surface tension and engaged binding sites per protein through a local order parameter, are provided in the Supporting Information document.