Liquid-Liquid Phase Separation of Patchy Particles Illuminates Diverse Effects of Regulatory Components on Protein Droplet Formation

Recently many cellular functions have been associated with membraneless organelles, or protein droplets, formed by liquid-liquid phase separation (LLPS). Proteins in these droplets often contain RNA-binding domains, but the effects of RNA on LLPS have been controversial. To gain better understanding on the roles of RNA and other macromolecular regulators, here we used Gibbs-ensemble simulations to determine phase diagrams of two-component patchy particles, as models for mixtures of proteins with regulatory components. Protein-like particles have four patches, with attraction strength εPP; regulatory particles experience mutual steric repulsion but have two attractive patches toward proteins, with the strength εPR tunable. At low εPR, the regulator, due to steric repulsion, preferentially partitions in the dispersed phase, thereby displacing the protein into the droplet phase and promoting LLPS. At moderate εPR, the regulator starts to partition and displace the protein in the droplet phase, but only to weaken bonding networks and thereby suppress LLPS. At εPR > εPP, the enhanced bonding ability of the regulator initially promotes LLPS, but at higher amounts, the resulting displacement of the protein suppresses LLPS. These results illustrate how RNA can have disparate effects on LLPS, thus able to perform diverse functions in different organelles.


Introduction
There is intense current interest in membraneless organelles, thought to be formed by liquid-liquid phase separation (LLPS) of protein-RNA mixtures. [1][2][3][4][5] Upon phase separation, these organelles emerge as droplets containing condensed component proteins, from the cytoplasm or nucleoplasm where the component proteins are dilute.
Membraneless organelles have been associated with many cellular functions. In particular, the condensation of component proteins may facilitate the assembly of complexes for biochemical processes, as exemplified by the nucleolus, an organelle with a primary function in ribosome pre-assembly. Component proteins usually contain intrinsically disordered regions (IDRs) and RNA-binding domains. The effects of RNA on LLPS have been investigated in many experimental studies, but conflicting results have been reported. 2,[6][7][8][9][10][11][12][13][14][15] The aim of this study was to use computation to elucidate potential roles of RNA and other regulatory components in LLPS and the underlying physical mechanisms.
A hallmark of LLPS (in contrast to, e.g., aggregation) is thermodynamic reversibility. That is, droplets can easily dissolve upon raising temperature or salt concentration and can reform when conditions are reverted. The reversibility comes because the two phases, droplet and dispersed, are in thermodynamic equilibrium. The two phases separate when the molecules can achieve the same low free energy by adopting two distinct types of configurations, with disparate concentrations and extents of intermolecular bonding. 16 The low free energy is achieved by distinct means in the two phases: e.g., high entropy (due to low concentration) in the dispersed phase but low enthalpy (from intermolecular bonding) in the droplet phase.
Raising temperature or salt concentration leads to reduced stabilization of the droplet phase by intermolecular bonding. Beyond a critical point, the distinction between the two phases vanishes, and hence droplets dissolve. The thermodynamic equilibrium between the two phases is characterized by a phase diagram, or binodal, which is specified by the concentrations in the two coexisting phases for conditions below the critical point. In general, the stronger the intermolecular bonding, the higher the 4 critical point.
Among the experimental studies on LLPS of protein-RNA mixtures, most reported that the addition of RNA promoted LLPS, indicated by either reduced threshold protein concentrations for droplet formation or increased critical salt concentrations. 2,7,8,[10][11][12][13] Reinforcing the promotional effects, both Saha et al. 12 and Smith et al. 13 found that molecules that competed for RNA binding with the droplet-forming proteins also suppressed phase separation. Intriguingly, Burke et al., 6 Zhang et al., 9 and Banerjee et al. 14  It should be recognized that the effects of RNA on protein LLPS are a reflection of intermolecular interactions rather than anything uniqueness about RNA. Indeed, similar promotional effects were observed in one case when RNA was replaced by heparin (which, like RNA, is an acidic polymer and thus presumably binds to the same sites on the droplet-forming protein) 10 and in another case, instead of RNA, the regulatory component was arginine-rich peptides (expected to bind to acidic tracks, rather than RNA-binding sites, of the droplet-forming protein). 11 LLPS also occurs for, and in fact, was first observed on, globular proteins. 17,18 The phase separation of droplet-forming proteins was suppressed, with decreased critical temperatures, when another protein was added. 19,20 For both disordered and globular proteins, LLPS was 5 promoted by crowding agents such as polyethylene glycol (PEG) or Ficoll. 7,8,21,22 Interestingly, in the latter study, bovine serum albumin (BSA) was found to suppress the phase separation of hnRNPA1 or its IDR fragment, but promoted, as if as a crowding agent, the phase separation of a protein-RNA mixture. Moreover, lysozyme exhibited a concentration-dependent dual effect on the LLPS of the IDR fragment of FUS (FUS IDR ). In short, RNA and other regulatory components can exert similarly diverse effects on the LLPS of disordered and globular proteins. This is not surprising, since the LLPS of disordered and globular proteins have a common physical basis, though their phase behaviors show characteristic differences. 23 Experimental studies have presented many observations on LLPS that beg for mechanistic answers. In principle, theoretical calculations and molecular simulations can provide these answers, but still face significant technical challenges. Polymer theories, 24,25 coarse-grained simulations, 26 and hybrids thereof 27 have been used to determine phase diagrams of chain molecules, in particular as models of intrinsically disordered proteins (IDPs), and to address several important questions. These include how phase boundaries are affected by charge patterns along the sequence, chain length, charge mutations, and salt. The study that relates most closely to the present one is by Lin et al., 25 who treated the phase separation of a mixture of two IDP sequences within a polymer theory. All these approaches have limitations and hence their predictions are necessarily qualitative. In particular, polymer theories often provide only a very crude account of excluded-volume effects. 23 Here we used Gibbs-ensemble simulations to investigate how regulatory components affect protein droplet formation, a pressing question that has received only scant attention in previous theoretical and simulation studies. While the promotion of LLPS by crowding agents can be understood from the perspective of excluded-volume effects, 28,29 when RNA or another regulatory component is attractive toward the droplet-forming component and hence is recruited into the droplet phase, the physical rules governing the foregoing disparate effects on LLPS are yet to be established. Toward that aim, we chose a relatively simple model for 6 protein-regulator mixtures, to calculate phase diagrams over a wide range of parameter values. To our knowledge, this is the first computational study devoted to uncovering the physical determinants underlying the roles of RNA and other regulatory components in LLPS. Our results qualitatively recapitulate the behaviors observed in the experimental studies, and provide a physical basis, in terms of the relative strengths of protein-protein and protein-regulator interactions, to reconcile the apparently conflicting effects of RNA on droplet formation.

Results
A simple model for protein-regulator mixtures. The droplet phase of protein-RNA mixtures, as true for LLPS in general, is stabilized by bonding networks of the molecules. 23 As alluded to above, we reason that the disparate effects on LLPS are generic instead of unique to RNA, and hence can be captured by models that possess common features of mixtures of proteins with RNA or other regulatory components, in particular, polyvalent interactions. Here we used patchy particles to model protein and regulatory molecules (Figure 1a). These particles, with attractive patches decorated on a hard core, which in our case is spherical, have been used to model proteins in LLPS. [30][31][32][33] The details of the protein and regulator models were chosen based on a number of observations regarding LLPS of protein-regulator mixtures. First of all, in many experimental studies where LLPS has been demonstrated, there exists a single "driver" protein that can form droplets on its own, whereas RNA or other regulatory components can modulate the phase boundaries. Correspondingly, interprotein interactions have to be strongly attractive (for droplet formation), whereas inter-regulator interactions are weakly attractive or even purely repulsive. The latter scenario is certainly true for RNA, due to the uniform anionic charges of the phosphates. Meanwhile protein-regulator interactions are expected to be highly dependent on the particular molecules involved. Lastly, all particles occupy volume so they experience steric repulsion, which, much like the attractive interactions 7 between patches, should play a role in determining phase diagrams, as demonstrated by effects of crowding agents.
In our models, protein-like ("P") and regulatory ("R") particles have 4 and 2 patches, respectively, to be denoted as ! ! , where S = P or R. The patches together cover the same fraction, !, of the surface areas of the hard cores, which have the same diameter ! for both types of particles. Each patch is a spherical disk, with the spanning polar angle ! ! given by Each patch is identified by the unit vector from the particle center to the patch center, denoted by ! with appropriate indices. The centers of the patches are positioned as a tetrahedron on a P particle but at the poles on an R particle. The greater spreading of the patches on a P particle gives it an advantage in forming bonds with multiple partner particles. A bond is formed when a patch on one particle is in contact with a patch on another particle.
The interaction energy between any two particles, i and j, has a directional square-well form, 34 where ! !" is the distance between the particle centers, ! !" is the unit vector between them, ! ! and ! ! denote the orientations of the particles. The square-well part has the usual form, where ! !" and ! are the strength and width, respectively, of the interparticle attraction. The interparticle attraction operates only when the patches of the two particles form a bond. The latter condition is specified by ! !" • ! !" > cos ! !" and −! !" • ! !" > cos ! !" , where α or β refers to any of the patches on particle i or j, and ! !" or ! !" is the corresponding spanning polar angle of the patch. Hence 8 Note that each patch can simultaneously form bonds with multiple partner particles.
The total energy of the two-component patchy particle system is the sum of all pairwise interactions between particles.
In the present work, ! was fixed at 0.7 and ! was fixed at 0.5!. We set ! !! = 0 (so the R particles experienced only steric repulsion between each other), and varied the strength of attraction between P and R particles, ! !" , over a wide range, from 0 to 1.5! !! .
Phase diagrams. For each ! !" , we used Gibbs-ensemble Monte Carlo simulations 35 to determine phase diagrams for the protein-regulator model system over a wide range of molar ratios between P and R particles. In these simulations, fixed numbers of P and R particles were placed in a box with a fixed volume and maintained at a constant temperature. The box was divided into two regions, and the volume and the particle numbers in each region were changed to achieve equality in pressure and chemical potentials. After reaching equilibrium, the two regions represent different phases, with distinct concentrations (denoted by ! with appropriate subscript and superscript) for the two types of particles ( Figure 1b). We refer to the low-and high-concentration phases by I and II. A binodal consists of coexistence concentrations, ! ! ! and ! ! !! , at a series of temperatures. For notational simplicity, hereafter all lengths and energies are, respectively, in units of the particle diameter ! and the strength of P-P attraction, ! !! (equivalent to defining ! = 1 and ! !! = 1). Likewise the temperature, T, is in We denote the ratio of total R and P particle numbers as X.
The effects of the R particles on the phase diagrams differ qualitatively depending on ! !" (Figure 2). They fall into three regimes, at low, moderate, and high ! !" values. The low-! !" regime is represented by the results at ! !" = 0 Here the effects of the R particles depend on the R to P molar ratio. With increasing X, T c initially increases, reaching a maximum at X = 0.35 (indicated by a parabolic gray arrow in Figure 2d). Then T c gradually decreases, going below the counterpart for the pure P system at approximately X = 1.3. The threshold P concentration for phase separation has a similar turnover behavior, reaching a minimum at X = 0.7 for T = 0.7 (Figure 2d inset). In contrast, ! ! !! is a monotonically decreasing function of X. When ! !" is further increased to 1.5, the maximal T c shifts to a higher X value of approximately 0.65 and promotion of LLPS (i.e., T c > the counterpart for the pure P system) persists over a wider range of X ( Supplementary Fig. S1c).
Partitioning of R particles in two phases. While Figure 2 displays the concentrations of the P particles in the two coexisting phases, the phase equilibrium is fully characterized by the concentrations of both types of particles. In Figure 3 we present the compositions of the two phases for the system with the four representative ! !" values at a given temperature. At each X, the compositions are shown as a pair of , connected by a black line (known as a tie line).
At ! !" = 0, the R concentrations in the droplet phase are miniscule for all X values, while those in the dispersed phase increase with increasing X (Figure 3a).
These R particles preferentially partition in the dispersed phase, as expected for volume-exclusion crowders. In doing so, they displace the P particles into the droplet phase, as indicated by higher and higher ! ! !! and lower and lower ! ! ! at increasing X.
The latter is already seen in Figure 2a as near symmetric broadening of the binodal.
Preferential partitioning in the dispersed phase occurs at ! !" up to 0.2, where steric repulsion of R particles still dominates.
With increasing ! !" , more and more R particles are recruited into the droplet phase. Whereas the R particles still prefer the dispersed phase at ! !" = 0.5, as indicated by tie lines with a negative slope (i.e., ! ! ! > ! ! !! ; Figure 3b), the preference switches to the droplet phase at ! !" = 1.0, where the tie lines now have a positive slope (Figure 3c). The preference of the R particles for the droplet phase further strengthens at ! !" = 1.35 ( Figure 3d). As more and more R particles are recruited into the droplet phase, more and more P particles are displaced from it. The latter trend is already seen in Figure 2c,d as significant leftward shifts by the high-concentration arm of the binodals. The simple reason for the displacement is that both P and R particles occupy volume so that the total concentration of particles that can be accommodated in the droplet phase is limited.
Correlation between T c and average number of bonds. The critical temperature for LLPS is determined by the strength of the bonding networks in the droplet phase. We measured the latter by the average number of bonds, n b , formed per particle at a given temperature. In Figure 4 we present T c and n b side by side, as functions of X. For all the four representative ! !" values, the dependences of T c on X qualitatively follow those of n b . Specifically, at ! !" = 0, both T c and n b increase with increasing X (though a turnover may be expected at very high X; Figure 4a), whereas at ! !" = 0.5 and 1.0 both T c and n b decrease monotonically with X (Figure 4b,c). At ! !" = 1.35, both T c and n b exhibit a turnover, with the maxima occurring at nearly the same X value (around 0.3).
Physical mechanisms for diverse effects of R particles. We are now in a position to elucidate the physical mechanisms behind the various effects of the R particles on phase separation, by examining how they affect n b . At low ! !" , as typified by the situation at ! !" = 0, the R particles preferentially partition into the dispersed phase, thereby displacing more and more of the P particles into the droplet phase. As ! ! !! increases, each P particle is surrounded by more and more other P particles, and hence can form more bonds. This explains why n b at T = 0.7 increases from 6.6 at X = 0 to 7.5 at X = 1.5 (Figure 4a). Therefore volume-exclusion crowders promote LLPS, as observed for both protein-RNA mixtures and globular proteins, 7,8,21,22 by preferentially partitioning into the dispersed phase and thereby passively displacing droplet-forming proteins into the droplet phase to form stronger bonding networks.
To better explain how the R particles affect n b at moderate and high ! !" , we further decomposed it. n b is the average of n b;P , the number of bonds formed per P particle, and n b;R , the counterpart per R particle, weighted by their populations (! ! !! and ! ! !! ): where ! !! = ! ! !! + ! ! !! is the total concentration of particles in the droplet phase.
Since a P particle can form bonds with both P and R partners (whereas an R particle can only have P partners), we further separated n b;P by the types of partner particles, In Figure 5 we present the decomposition of n b as well as ! ! !! /! !! at ! !" = 0.5,

12
bonding ability than the P particles (recall that ! !! = 1 by definition). Indeed, when surrounded by pure P particles, a P on average forms 7 bonds [! !;! at X = 0 in Figure   5a, at T = 0.68), while an R forms only 4.5 bonds (extrapolation of ! !;! to X = 0). As X increases, more and more R particles are recruited into the droplet phase, and hence they weight down the overall n b . In addition, at the highest X value (i.e., 0.43) studied, ! !;! is significantly lower (to 6.5 bonds) than the counterpart in the pure P system.
The latter comes about because, due to the displacement of P particles by R particles, the loss in ! !;!!! , the number of bonds with P partners, is not fully compensated by the gain in ! !;!!! , the number of bonds with R partners.
The behavior of n b at ! !" = 1.0 is similar to that at ! !" = 0.5, but shows subtle differences ( Figure 5b). The bonding ability of the R particles is now higher. Still, when surrounded by pure P particles, an R still lags behind a P in the number of bonds formed (6.5 for R versus 7 for P), due to less spreading of the patches on the R particles. Moreover, ! !;! now decreases significantly with increasing X, due to decreasing ! ! !! and hence a less number of bonding partners. The resulting widening gap between ! !;! and ! !;! , coupled with the increasing R population, accounts for most of the decrease in n b with increasing X. Lastly, ! !;! itself decreases with increasing X, again because the gain in ! !;!!! does not keep pace with the loss in ! !;!!! . In short, at moderate ε RP , the R particles, by partitioning into the droplet phase and displacing P particles, actively weaken bonding networks there and thereby suppress LLPS.
At ! !" = 1.35, the stronger R-P attraction (relative to the P-P attraction) enables an R to form more bonds than a P when each is surrounded by pure P particles (8 for R versus 6.5 for P at T = 0.7; Figure 5c). With the enhanced bonding ability for R, as X increases from 0 and more R particles displace P particles, the gain in bonds for a P formed with R particles is greater than the loss in bonds with other P particles ( Figure   5d, left). Consequently, ! !;! , the average number of bonds formed by a P, initially increases with increasing X. After X = 0.25, ! !;! levels off. Meanwhile ! !;! decreases monotonically with X due to the displacement of partners (i.e., P particles) 13 by R particles (Figure 5d, right). At around X = 0.3, the weighted average of ! !;! and ! !;! , i.e., n b , thus reaches maximum, and hence here the R particles have optimal ability in promoting LLPS. As X increases further and further and more and more P particles are displaced from the droplet phase, the R particles cannot form enough bonds so eventually they suppress phase separation.

Discussion
We have calculated the phase diagrams of a model protein To elucidate the underlying physical mechanisms of the effects of regulatory components, we have found it useful to characterize bonding networks in the droplet phase, specifically, by the average number of bonds per particle and its decomposition by species. In general, molecular interaction networks, or transient clusters, are crucial to the physical understanding of dense macromolecular systems, such as those in cellular environments. 16,36 Our results qualitatively recapitulate the behaviors observed in experimental studies of protein-crowder, protein-protein, and protein-RNA mixtures, and provide a physical basis to reconcile conflicting reported results for the effects of RNA on protein droplet formation. Crowding agents such as PEG and Ficoll are found to promote LLPS for protein-RNA mixtures and globular proteins. 7,8,21,22 Our calculation results with low ! !" show that the promotional effect of such crowders comes because they preferentially partition in the dispersed phase, thereby passively displacing droplet-forming proteins into the droplet phase to strengthen the bonding networks there.
The critical temperature for LLPS of γD crystallin decreased upon adding βB1 crystallin. 19 Similarly, droplet formation of a monoclonal antibody and hnRNPA1, respectively, was suppressed by human serum albumin 20 and BSA. 22 The changes in T c and the tie lines of the protein mixtures in the first two studies are similar to our results obtained at moderate ! !" . In this case, the regulatory protein has weaker attraction toward the droplet-forming protein than the latter toward itself. Still, the regulatory protein is recruited to the droplet phase, and thereby displaces the droplet-forming protein and weakens the bonding networks there. Both the recruitment and displacement features identified by our study were directly visualized by fluorescence microscopy for the hnRNPA1-BSA systems, whereas PEG and Ficoll, also as expected, were excluded from the droplet phase. 22 Displacement occurs simple because there is a limit to the total concentration of macromolecules in the droplet phase.
Promotional effects of RNA on LLPS have been reported in a number of studies. 2,7,8,[10][11][12][13] We explain these results by suggesting that the RNA molecules in these studies have stronger attraction toward the droplet-forming proteins than the latter toward themselves, and therefore actively strengthen bonding networks in the droplet phase. This suggestion is quite reasonable, since the droplet-forming proteins all contain specific RNA-binding domains to achieve strong bonding with RNA. We further hypothesize that the promotional effects in these studies were obtained when the amounts of RNA added were sufficiently low, and predict that adding higher amounts of RNA would eventually lead to LLPS suppression. Indeed, such dual effects have been observed in three studies. 6,9,14 The present study reveals that the suppression at high RNA-to-protein molar ratios arises from the displacement of proteins by RNA from the droplet phase, leading to an insufficient level of proteins to maintain the strength of the bonding networks of the protein and RNA molecules.
This mechanism is distinct from those by Zhang et al. 9 and Banerjee et al. 14 , who speculated that RNA would either screen interprotein electrostatic interactions like diffuse salt ions or invert the protein net charge by strong binding to the protein surface. While these specific effects of RNA may exist, our mechanism is generally applicable to RNA and any other macromolecular regulators.
As support for our contention that the findings in the present study are applicable to both rigid and flexible molecules, a concentration-dependent dual effect was also found for lysozyme, a globular protein, on the LLPS of FUS IDR . 22 In comparison, BSA did not show any sign of promoting the FUS IDR LLPS. We can now explain the different effects of lysozyme and BSA in terms of their relative strengths of attraction toward FUS IDR . That is, we predict that FUS IDR is more strongly attracted to lysozyme than to BSA. Interestingly, a recent small-angle neutron scattering study indicated exactly such differential strengths of attraction of the latter two proteins toward FlgM, an intrinsically disordered protein. 37 While two regulatory proteins (lysozyme and BSA) can have divergent effects (active promoter vs. active suppressor) on the phase separation of a given protein (FUS IDR ), a given regulatory protein can also have distinct effects on the phase separation of two different proteins. The latter was observed with BSA acting as an active suppressor for the LLPS of hnRNPA1 but as a volume-exclusion crowder for the LLPS of a protein-RNA mixture. 22 Again, the distinct effects can be attributed to the differential strengths of attraction of the regulatory protein toward the droplet-forming proteins. We note that the strengths of intermolecular interactions can be directly measured by the second virial coefficient or the binding affinity, determined at concentrations below the threshold for LLPS. 13,15 In another intriguing report, 15  Our study has identified the relative strengths of protein-protein and be induced by adding RNA. We hope that these and other predicted trends in Figure 6 will motivate experimental tests.

Methods
The phase diagrams of the protein-regulator model system for ! !" (strength of P-R attraction) between 0 and 1.5 and X (R to P molar ratio) between 0 and 1.5 were determined from Gibbs-ensemble Monte Carlo simulations. 35 In most simulations, the total number of particles was 512 and the total concentration, ! ! , was 0.3, a setup where the binodal for the pure P system has been obtained previously; 34 our simulations at X = 0 reproduced this result. Some simulations were also carried out with a total of 1000 particles (while preserving ! ! ); the resulting coexistence concentrations were unchanged. In other simulations ! ! was changed to 0.1 and the results are qualitatively consistent with those at ! ! = 0.3. Each binodal was obtained from simulations at a series of temperature ranging from 0.65 to 0.80. At each temperature, the first 1 to 1.5 million cycles were discarded to ensure equilibration, and the subsequent 2 to 3 million cycles were used for data collection. Each cycle consisted of 500 attempts of particle displacement, 500 attempts of particle exchange, and 5 attempts of volume exchange.
The critical point for each binodal was determined by fitting the temperature dependence of the coexisting P concentrations, ! ! ! and ! ! !! , to the law of rectilinear where ! ! ! is the critical protein concentration, and A and B are other fitting parameters. Each fit was performed first by including data at all available temperatures and then gradually removing data that were closer to T c , until the fitted values of ! ! ! and T c were stable.
Data availability. The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.     are shown as full-size spheres in red and blue, respectively, with bonds in gray;

Figure Captions
nonbonded P and R neighbors (2 Ps in the left and 4 Rs in the right) are shown as full-size spheres in cyan and magenta, respectively. The latter coloring scheme is also used for the more distant P and R particles, shown with radius reduced to 40%. Figure 6. Expected effects of regulatory ("R") particles on the droplet formation of protein ("P") particles. Suppression of LLPS occurs at moderate P-R attraction (blue color), whereas promotion of LLPS occurs at either very weak or strong P-R attraction (red color). The promotion can turn into suppression of LLPS at high R to P ratios.