Generating lattice representations of multivalent proteins

For a given linear or branched multivalent protein, we first choose a suitable mapping between the protein degrees of freedom and a lattice representation. The conformational space is a simple cubic lattice with periodic boundary conditions used to mimic a macroscopic system. Phase transitions represent the collective effects of large numbers of molecules, and simulations have to include at least 10³–10⁴ protein molecules to observe facsimiles of these collective transitions in finite sized systems [82]. Further, we need to be able to test for the effects of finite size artefacts and this requires a titration of the effects of varying the numbers of molecules. Accordingly, the lattice has to be large enough to accommodate at least 10³ molecules of each type for the most dilute concentrations. Often, we might need to increase the number of molecules to be of O(10⁴). Accordingly, a one-to-one mapping between the protein degrees of freedom and a lattice representation would lead to a computationally intractable model. Instead, we adopt system-specific coarse-graining approaches, whereby the coarse-graining is guided by a priori rigorous or phenomenological knowledge of the identities of stickers versus spacers. For disordered proteins, the stickers within disordered regions often correspond to single amino acid residues or short linear motifs. For multivalent folded proteins, the stickers are either an entire protein domain or sectors on domain surfaces [28, 29]. Residues corresponding to spacers may either be modeled explicitly, where one or more spacer residues are modeled by a single bead on the lattice site, or be modeled as phantom tethers, where the intrinsic lengths of tethers are calibrated in terms of the numbers of lattice sites [28, 29]. In both cases, the tethers can stretch, bend, and rotate and these degrees of freedom contribute to density inhomogeneities that are the result of altered patterns of inter-sticker interactions.

LASSI and bond fluctuation models

The structure of LASSI is inspired by the bond fluctuation model (BFM) for lattice polymers [83]. This is a general lattice model for simulations designed to extract equilibrium conformational distributions and dynamical attributes of polymers in dilute solutions as well as dense melts. There are two versions of the BFM viz., the Carmesin-Kremer BFM or CK-BFM [84] and the Shaffer BFM or S-BFM [83]. Both models are based on the use of simple cubic lattices, which discretizes the conformational space for polymers.

In the CK-BFM [84], each repeating unit or monomer within a polymer is modeled as a 3-dimensional cube where the 8 corners of the cube occupy lattice sites and bond vectors connect pairs of monomers. Overlap of monomers is associated with an energetic penalty, and each bond vector can have up to 108 distinct directions. The choice of bond vector set encodes the geometry of the polymer and places constraints on the bond lengths and bond angles. All other interactions are governed by the inter-monomer potentials, and evolution of the system through conformational space is driven by changes to the overall potential energy. In contrast, the S-BFM places each monomer on a single lattice site. Covalently bonded monomers are connected by bonds that are constrained to be of three types, leading to chains that have bonds of length 1, or in units of lattice size. Monte Carlo moves with suitable acceptance criteria can be designed for both types of BFMs. The simulations are used to generate equilibrium conformational distributions of lattice polymers in either dilute or dense phases. The move sets control the overall polymer dynamics and the acceptance of different types of moves and the calculation of correlation functions allows one to compute dynamical quantities for lattice polymers [83]. If we were to use either of the established BFMs without modification, then each amino acid residue would be modeled as a monomer, and such an approach would be useful when the identities of stickers and spacers remain ambiguous. This approach underlies a different simulation engine known as PIMMS [43].

LASSI is a generalization of the S-BFM that also adapts features of the CK-BFM. Given a choice of the mapping for coarse-graining, each multivalent protein is described as a chain of non-overlapping monomers viz., beads that occupy sites on a 3-dimensional cubic lattice. Note that the choice of a single site per bead is similar to that of the S-BFM, although the bead, which is a sticker or spacer monomer, need not be the monomeric unit, i.e., an amino acid residue in the case of proteins. Each sticker monomer is linked to its adjacent sticker on the chain via either a phantom tether or a set of spacer beads that occupy individual lattice sites [28, 29]. A spacer / tether length of unity implies that adjacent monomers are within lattice units of one another (Fig 1D). The choice of the spacer length will be sequence-specific or more precisely, specific to the architecture of the protein of interest.

Inter-monomer (sticker-sticker, sticker-spacer, and spacer-spacer) interactions are modeled as contact-based pairwise interactions. A sticker monomer can bind to another sticker monomer that occupies an adjacent lattice site with an interaction energy that depends on the types of both monomers. Monomers are considered to be adjacent to one another if they are within a lattice distance of . By this criterion, each lattice site occupied by a sticker monomer will have 26 adjacent lattice sites. This is reminiscent of the interaction geometry of a CK-BFM for each monomer. In the current implementation of LASSI, the interactions are mutually exclusive, implying that a sticker cannot interact simultaneously with more than one other sticker, even though there are 26 adjacent sites that the interaction partner can occupy. If the sticker in question is already engaged in another inter-monomer interaction with stickers or spacers, then the unoccupied sites of the sticker will be unavailable for interaction. The combination of the geometry of the interaction sites per monomer and the single occupancy constraint leads to anisotropic interactions between sticker interactions. This feature is unique to LASSI and it is not incorporated in other variants of BFMs; this allows us to deploy LASSI for modeling heteropolymeric systems. In the context of LASSI, we note that stickers are distinguished by their ability to participate in anisotropic or isotropic interactions. In contrast, explicitly modeled spacer sites only participate in isotropic interactions with other spacer or sticker sites. Furthermore, the interaction strengths involving spacers are typically weaker than those involving stickers. However, it is worth emphasizing that these distinctions only matter inasmuch as LASSI allows us to capture a numerical instantiation of the stickers-and-spacers model. For simplicity, one might simply think of LASSI as a model that has sites that are differentiated by whether or not they can involve themselves in anisotropic interactions, by their intrinsic site valence (a variable that we do not titrate in this work), and by the comparative magnitudes of site-site interaction strengths.

Setup of simulations

A system with n multivalent proteins is in reality an n+1 component system since the solvent is the implicit component. In LASSI, sites that are not occupied by protein units automatically represent solvent sites. Although the interaction potentials do not explicitly include terms between solvent and protein sites, the effective interaction strengths between pairs of protein units represent an averaging over protein-protein, protein-solvent, and solvent-solvent interactions. The solvent sites, i.e., the sites that are not occupied by protein units, represent contributions from the solvent to the overall translational and mixing entropies. Simulations are initiated by randomizing the positions of protein units, subject to the constraints of chain connectivity.

The parameters that are set at the start of each LASSI simulation include the total number of molecules n_i of type i and the size of the lattice L, from which we can calculate the total number n of all protein components and the concentration or number density of each protein c_i = n_i/L³. The setup also includes stipulations for the architectures of each protein such as specification of the number of monomers per chain, the overall topology of each protein (linear vs. branched), the lengths of spacers, and the types of spacers (implicit / phantom vs. explicit) [28–30, 79, 80]. The number of monomers per molecule will equal the sum of the number of stickers and spacers if spacer residues are modeled explicitly. Alternatively, if spacers are modeled as phantom tethers, then the number of explicitly modeled monomers will equal the number of stickers. Specification of the energetics of the system includes specification of the simulation temperature in normalized units, homotypic and heterotypic interaction strengths between pairs of stickers, the energetic cost for the overlap of stickers, and the interaction strengths between sticker and spacer sites if the spacers are modeled explicitly.

Design of monte carlo move sets

Our goal is to compute architecture-specific phase diagrams for systems comprising of one or more types of linear or branched multivalent proteins. This requires a simulation strategy that enables the sampling of the full spectrum of coexisting densities and networked states for multivalent proteins. Accordingly, the conformations of randomly initialized systems of proteins on a simple cubic lattice are sampled via a series of Markov Chain Monte Carlo (MCMC) moves that are designed to ensure efficient sampling of changes in protein density and networking while maintaining microscopic reversibility. We have developed and deployed a collection of moves and these are described below.

Monte carlo sampling with biases

In LASSI, we have independent contributions from two main energetic sources. Monomer units are not allowed to overlap, and this can be described by a position-dependent energy E_pos where E_pos = 0 or ∞. On the other hand, inter-monomer pairwise interactions also contribute to the total energy, and E_rot denotes the sum over all of the effective pairwise inter-monomer interaction energies. The subscript “rot”(rotational) indicates the fact that for a pair of nearest neighbor stickers their interaction energies are actually governed by their mutual orientations. Accordingly, the total system energy in a specific configuration i is written as: (1) The equilibrium probability associated with configuration i is given by the Boltzmann distribution as: (2) In Eq (2), β is the inverse of the simulation temperature in units of the Boltzmann constant (effectively, k_B = 1 energy unit / temperature unit). The frequency with which a transition from configuration i to j is proposed will be governed by the elements T_ij of the targeting matrix T. The proposed transition is accepted / rejected based on the elements A_ij of the acceptance ratio matrix A. A MCMC move that transitions the system from configuration i to j defines a flow in configuration space and this flow has to satisfy microscopic reversibility: (3) If the targeting matrix is symmetric, then T_ij = T_ji and the acceptance ratios such as those prescribed by Metropolis et al. [85] will ensure the preservation of microscopic reversibility. However, if biases are incorporated into the targeting matrix, which is often essential to enhance the sampling of configurations that contribute to density inhomogeneities and the making / breaking of bonds in dense networks, then the elements of the acceptance ratio matrix have to be designed to ensure the preservation of microscopic reversibility. We deploy a general strategy of using biased moves to enhance the sampling of different mutual orientations among pairs of stickers. The incorporation of these orientational biases is accounted for by modifying the acceptance criterion of Metropolis et al. [85] whereby each element of A is written as: (4) For a symmetric targeting matrix, we recover the standard acceptance ratio of Metropolis et al. [85] viz., (5) Since the moves within LASSI generally involve orientational biases, the elements T_ij are rewritten in terms of a Rosenbluth weighting factor W_j [86, 87] whereby: (6) Substituting (6) into (4) leads to: (7) The specific form for the weighting factors W_i will depend on the type of move because the extent of asymmetry in the targeting matrix will depend on the nature of the bias incorporated into the biasing move that proposes a transition from i to j. The specific forms for weighting factors are discussed in the context of the move types that are introduced next.

Rotational moves

A monomer is in an associated or a dissociated state and in the associated state it has a specified binding partner. This defines the rotational state of a monomer. To change the rotational state, we randomly pick a monomer from the system, and exhaustively sample all 26 adjacent sites to construct a list of potential binding partners. The rotational state of the monomer is changed, at random, by drawing a random integer k from a uniform distribution between [0, b], where b is the number possible binding partners available to the monomer. If k = 0, the monomer is set to be in a dissociated state. Otherwise, the k^th candidate bond (reversible physical crosslink) is formed and the state of the monomer is set to be in an associated state. If the monomer cannot be involved in a rotational interaction, as would be the case for an explicitly modeled spacer, the rotational move is rejected outright. The accessible volume for rotational interactions is within a cube of unit volume centered on the randomly chosen monomer (see S1 Fig for a 2-dimensional representation), and hence, each sticker monomer will have at most b_max = 3³–1 possible sites as neighbors. Since this number is not large, we sample all 26 possible interaction sites.

Local moves

This move serves as the basic unit of local displacement of monomers–be they stickers or spacers. A randomly chosen monomer is moved from position r_i to r_j = r_i + Δr. Acceptance of the move is predicated on the move not leading to an overlap with a site occupied by another monomer and the satisfaction of linker constraints. The choice for Δr is made by uniformly sampling each component from the interval [–2,2] such that (shown in S2 Fig) moves if the selected monomer is in the interior of a molecule [88], and they become analogous to end rotation moves if end monomers are selected [89].

Local moves have a rotational bias in LASSI and the Rosenbluth factor is calculated as follows. Starting with Eq (7), we shall designate the chosen monomer by index k. In configuration i, assume that monomer k has a binding partner of index l. Typically in coarse-grained systems there are a finite number of unique monomer types, and thus it is more efficient to simply define interaction energies between different monomer types than between all monomer pairs. The energy associated with the bond between monomers k and l is written as ε_t(k)t(l), where t(x) indicates the type of monomer x. The local move causes a change in binding partner, whereby the monomer k now binds to monomer m. The local move leads to a bond swap that causes a change in rotational energy, which is written as: (8)

Use of Eq (8) in Eq (6) leads to: (9)

In Eq (9), each Rosenbluth weight has an additional index in the subscript to indicate that the change in configuration is achieved by a change in the binding partner for the monomer k. To accelerate the creation of density inhomogeneities in supersaturated systems and facilitate the making and breaking of networks, we decompose W_i;k as: (10)

The two terms in Eq (10) respectively represent the contributions to the Rosenbluth weights for monomers in associated (a) and dissociated states (d). First, we calculate the weight factors for the interacting monomers as a partition function over all nearest neighbor contributions, such that: (11)

In Eq (11), the summation runs over all potential binding partners l (nearest neighbors) for the monomer k. To illustrate how the Rosenbluth factor is calculated, we assume that the system has only one type of interaction with the pairwise energy designated as ε. If the number of nearest neighbors for monomer k in configuration i is designated as N_k;i, then the Rosenbluth weight factor in Eq (11) becomes: (12)

Setting to incorporate a bias towards associated states, and using the simplification that leads to Eq (12), we rewrite Eq (9) as a definition of the acceptance criterion for the move from configuration i to j via local move involving monomer k as: (13)

This choice for the acceptance criterion ensures that detailed balance is preserved while enhancing the sampling of configurations characterized by the breaking of old bonds and the forming of new ones.

Reptation–or slithering snake–moves

In dense configurations, it becomes difficult to realize large-scale translational or rotational motions of polymers. The slithering snake move is a Monte Carlo instantiation of reptation as first conceived by de Gennes [90]. In this move, a chain is chosen at random, and the monomer at one end of the target chain is moved to a new position. The remaining monomers within the target chain are then successively moved such that monomer m along the chain moves into the previous position of monomer m–1 (S3 Fig). This move relies on an inherent symmetry of chain molecules, because bond lengths between monomers are the same; if one swaps monomers across chains, the identity of the chain remains invariant. However, this move cannot be used if the molecule has heterogeneous bond lengths or if it is a branched polymer.

The reptation move is rotationally biased, and this is true for every monomer in a chain. The bias is independent for each monomer and accordingly, the Rosenbluth factor for a single reptation move can be calculated from the Rosenbluth factors for each monomer-specific local move. In configuration i we obtain: (14)

In Eq (14), the product runs over all monomers m within the chain of interest. The acceptance criterion for a reptation move takes the form: (15) The inclusion of the bias for every interacting monomer, rather than just the end monomers, is to emulate how a real transiently bonded polymer would slither along its contours. Note that for strict detailed balance, the Rosenbluth factors for the two end monomers should be calculated and added to the acceptance criterion, but the current implementation of LASSI uses the first trial position that satisfies the position constraints.

Double pivot moves

These moves swap a part of a chain with the corresponding part of another chain of the same type. A monomer is picked at random; it is denoted as i_m, where i is the monomer index within a chain, and m is the chain index. A search is then performed around the monomer within a prescribed distance for a monomer within the same type of chain. The requirement for the search is that the monomer of interest be one index ahead along its own chain, (i+1)_n. Next a check is performed to ensure that the distance between i_m and (i+1)_n is within the bond constraint connecting i_m to (i+1)_m, and that the distance between (i+1)_m and i_n is within the bond constraint for i_m and (i+1)_m. Each monomer (i+1)_n that satisfies each of these constraints is stored and one of these is randomly picked for the double pivot move (S4 Fig). The move is always accepted if there is a candidate because only connectivity changes unless bonds are modeled using elastic springs.

The purpose of this move is to engender large configurational changes in dense polymer melts, which approximate the dense phases formed upon phase separation. In dense regions the rate of acceptance of local moves decrease precipitously. At high enough densities, polymers become entangled and local moves reduce to slithering-snake moves and polymers are restricted to motions along tubes around one another [91, 92]. Therefore, rather than physically moving polymers to create a change in configurations, we incorporate move sets that break and make bonds while ensuring that monomers do not overlap, and that bond constraints are always satisfied. If two chains are close enough to each other that the bonds between two monomers can be swapped, then such the double pivot move results in a large configurational change for both chains, and for the system.

Chain and cluster translation moves

The chain translation move is designed to move single chains while forming new bonds at the proposed location. This move attempts to translate the center of mass of a chain i from r_i to r_j = r_i + Δr where and L is the size of the simulation cell. Multiple trial displacements are proposed until a trial position that does not result in steric clashes for the entire molecule is generated. The move is then attempted. As with the slithering snake move, each monomer in the molecule that is translated will have an orientational bias. Accordingly, the Rosenbluth factors are calculated as in Eqs (14) and (15). The translational move results in large displacements for single chains and correctly biases the system for efficient sampling of configurations with alternate interaction patterns.

Translational moves can also be applied to clusters of molecules. A connected cluster refers to a collection of unique chains connected via rotational interactions. A proposed move only results in a translation, and the move is readily accepted if there is no steric clash. Since no new physical bonds are created at the proposed location, the cluster remains invariant and the move is accepted. Naively this move might seem unnecessary as this move simply moves clusters around. However, once a physical bond has formed between two molecules, it is highly unlikely for any of the non-cluster translation moves to move the centers-of-masses of clusters closer together.

Considerations for setting move set frequencies

The structure of each move set serves as a guide for selecting an optimal set of frequencies. This leads to a set of heuristics that are as follows: (i) in the cluster move we pick a random chain from the system, perform a networking analysis on that chain, and then propose a displacement of the cluster. As the cluster size grows it is more likely that a randomly picked chain will be part of the largest cluster which itself will result in a steric clash after the proposed move. Therefore, the frequency of the cluster move should be low, if not the lowest, in the entire set. (ii) In the translational move, we pick a random chain from the system for translation; as the size of the largest cluster increases it becomes less likely for a proposed translation move to be accepted. However, unlike the cluster move the translation move is rotationally biased and thus results in new interactions being formed. Hence, translational moves enable single-molecule to cluster-surface interactions. Therefore, this move should be proposed more often than the cluster move, although not as often as rotational or local moves. (iii) The rotation move is computationally inexpensive and it enables the switching of physical bonds and should thus be proposed fairly frequently. (iv) Similarly, local moves and slithering snake moves are also rotationally biased, and they help with the local rearrangements of physical bonds. Local moves are the primary route to enable local conformational changes, and to enable local physical bond rearrangements. Therefore, local moves should be proposed most frequently. The slithering snake move is particularly effective because it allows for large local physical bond rearrangements in dense configurations. Thus, this move should also be proposed frequently, less so than local moves but more so than translation moves. Note that in a system where some molecules are non-linear or have heterogeneous linker lengths, the frequency would need to be higher since the move is rejected if an incorrect molecule is picked at random. (v) The double pivot move allows for large-scale changes to conformations within dense configurations and accordingly, this move should be proposed more frequently than both cluster and translation moves. One can track the acceptance ratios of each move over a very rough initial sweep across the relevant system parameters. Moves that are always rejected do not enable any changes in configuration and only add computational costs. Therefore, the frequency for that particular move should be lowered. This is especially the case for the cluster move in high-density systems.

Identifying phase boundaries using measures for density inhomogeneities

In order to detect the onset of phase separation, we can calculate excess chemical potentials using the Widom particle insertion method [93] and equalize these chemical potentials across distinct phases. This process requires a priori knowledge of the densities of both phases. An efficient variant of this approach, based on fast Fourier transforms, was recently developed and deployed by Qin and Zhou [61]. They demonstrated their method for calculations of liquid-liquid coexistence curves for a patchy colloid model of γII-cyrstallin. Given that LASSI simulations are lattice-based, we instead rely on properties of pair distribution functions that help us diagnose the onset of phase separation and compute phase boundaries. Pair distribution functions are helpful because phase separation is the result of the system partitioning into phases of different densities. The pair distribution, which is a reduced-dimension partition function, serves as a rigorous thermodynamic and structural measure of the average local density and inhomogeneities of density. To first order, the density fluctuations are quantified by averaging over all orientations. Accordingly, the pair distribution function can be converted to a radial distribution function that allows us to probe local densities and local structural organization of molecules around one another. However, normalization of the pair distribution function requires some caution. The system contains polymer molecules and using a prior distribution that assumes an ideal gas of the chain monomers to normalize the pair distribution function is problematic because it does not accurately capture the effects of non-idealities due to chain connectivity. We leverage the efficient sampling of polymer fluids in LASSI and obtain suitable prior distributions by simulating the system of interest in the absence of sticker-sticker interactions.

The pair distribution function P⁽²⁾(r) quantifies the equilibrium distribution of distances between chain monomers, where r is the inter-monomer distance. If denotes the prior pair distribution function calculated from simulations where the inter-sticker interactions are ignored (see Fig 3A), then the normalized radial distribution function is written as: (16)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Distribution functions used for calculation of density inhomogeneity.

The data shown are obtained from 5 independent simulations for the A_n-B_n system with total protein concentration c = 6.89×10⁻⁵ voxels^-1 and reduced temperature T* = 0.383. Error bars indicate standard deviations. (a) Pair distribution functions P⁽²⁾(r) and P₀⁽²⁾(r), where the former is from the interacting system and the latter from the non-interacting system with chain connectivity (prior pair distribution function). Note that P⁽²⁾(r) shows two peaks, the first of which indicates dense phase formation. (b) Radial distribution function . This captures the droplet formation by a sharp and broad peak in the beginning. The inset shows r_b where intersects the line corresponding to line, delineating between the dense and solution phases. The global density inhomogeneity measure, , is obtained by integration of absolute deviation of from 1.

https://doi.org/10.1371/journal.pcbi.1007028.g003

The function is a direct measure of the local density of the protein of interest. Since LASSI uses periodic boundary conditions, the maximal inter-monomer distance is . Given this normalized , we note that if the system has short-range ordering as in a canonical liquid, the radial distribution function will oscillate around unity but eventually approach one as r → ∞. Conversely, if the system undergoes a density transition, will have two distinct spatial regimes (Fig 3B): for small r, will be characterized by a tall and broad peak such that until intersects the line; this region corresponds to the dense phase and we shall denote the value of r at this intersection to be r = r_b. For r > r_b, will be between 0 and 1, and for lattices that are large enough to avoid finite size artefacts, will converge to a value lower than one and this corresponds to the density in the dilute phase region. Furthermore, can be used to estimate the densities within the dense and dilute phases.

To quantify the global density inhomogeneity we introduce a simple measure, , which is calculated as follows: (17) In Eq (17), V(x), the volume element for the normalized radial distance x = r/L defined as: (18) (19) If ≈ 0, the global density inhomogeneity in the system is small and this will be characteristic of a single homogeneous phase dominating the simulation volume. As increases beyond zero, the system accommodates density inhomogeneities. We construct coexistence curves using a cutoff value of = 0.025, which is universal to all systems, to delineate between a homogeneous system and one that has undergone phase separation.

Quantitative assessments of finite size effects

The pair distribution function is central to our calculation of density inhomogeneities and constructing coexistence curves for a system of multivalent proteins simulated using LASSI. At the start of this section we emphasized the importance of including 10³–10⁴ distinct molecules within the simulation cell in order to avoid finite size artefacts. Prior to presenting detailed results that mimic specific systems, we present an analysis of finite size effects that we will confront if the requisite numbers of molecules are not included in the simulations. The data we present are for simulations of mimics of the protein FUS, specifically the A_n + B_n system introduced in the results section. The phenomenological mapping of this protein architecture onto a cubic lattice is discussed at the start of the results section. Here, we present an analysis that makes a crucial technical point about finite size effects.

First, we start with simulations for ideal polymers. The data shown in Fig 4 plots the pair distribution function extracted for simulations of ideal models of FUS-like proteins. Results are shown for simulations that use 20 chain molecules of A_n + B_n as an example of a small system. These results are compared to those from simulations with 100, 200, 1000, 1500, 2000, 3000, and 4000 A_n + B_n molecules, respectively. The pair distribution functions have a self-similar character and this is revealed by plotting for all of the simulations, where r* is the reduced distance that accounts for the fact that for a similar concentration, the simulation cells are made larger (higher values of L, which is box size) as the numbers of molecules increase. This analysis shows that even for a truly ideal system, the smallest simulation comprising of only 20 molecules will generate noisy estimates of the pair distribution function. This clearly demonstrates the problems inherent to small systems where finite size effects are accentuated. Interestingly, for the ideal chain system, all simulations with 100 or more molecules yield similar pair distribution functions as assessed in Fig 4B.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Assessment of finite size effects in simulations of ideal, non-interacting chains.

(a) Pair distribution functions computed in terms of the spatial separation between chain units. The distributions are maximal at r = L/2, where L is the size of the simulation cell for a given system. Note that L increases as the number of polymers in the system increases. With the exception of the smaller systems, the ideal chains show self-similar behavior for different system sizes. (b) The data plotted in panel (a) are re-plotted in terms of the scaled variable where L_i is the size of the simulation cell for boxes with i molecules.

https://doi.org/10.1371/journal.pcbi.1007028.g004

Next we assessed the impact of finite size effects with all of the terms in the potential being included in the simulation. There are three columns, one each for different values of the reduced temperature T*, in Fig 5. As discussed in the results section, these values of T* place the system of interest in the two-phase regime, with the quench depth into the two-phase regime increasing as T* increases. The first row (a to c) of Fig 5 shows the same data as Fig 4A while the second row (d to f) in Fig 5 show the normalized data like Fig 4B. Each panel shows eight unnormalized pair distribution functions, one each for the systems with 20, 100, 200, 1000, 1500, 2000, 3000, and 4000 molecules, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Assessment of finite size effects in simulations with real interacting chains.

Panels (a), (b), and (c), respectively are the real chain equivalents of panel (a) in Fig 4 computed for three different simulation temperatures that represent three different quench depths of the system into its two-phase regime. Panels (d), (e), and (f), respectively are the real chain equivalents of panel (b) in Fig 9 computed for three different simulation temperatures that represent three different quench depths of the system into its two-phase regime.

https://doi.org/10.1371/journal.pcbi.1007028.g005

The onset of phase separation should be manifest by the presence of a trough located between two peaks in the profiles for P⁽²⁾(r*). This is evident for T* = 0.267 (Fig 5, panel (c)) for all systems providing the numbers of molecules are greater than or equal to 10³. This qualitative trend is preserved even for T* = 0.217, although sampling difficulties in large systems become obvious in the noisy estimates for P⁽²⁾(r*). At the reduced temperature of T* = 0.167 we confront two problems: The small systems where the numbers of molecules are less than 10³ cannot support the distinction between a proper dense phase that coexists with a dilute phase. This behavior is similar to that observed for lower quench depths i.e., higher simulation temperatures as shown in panels (b) and (c) of Fig 5. However, as the system size grows, an additional problem arises and this has to do with large clusters becoming frozen, and thus inhibiting the achievement of equilibrium. This is evident from the pair distribution functions shown in panels (a) and (d) of Fig 5 for systems where the numbers of molecules exceed 10³. To overcome this broken ergodicity and obtain reliable converged pair distribution functions, we need additional biasing potentials and temperature sweep approaches used recently [43] to break up frozen clusters and enable their coalescence. In the results that we report here, we use the system size titration to identify the reduced temperatures below which broken ergodicity becomes evident. We do not include data from these simulations in our constructions of phase diagrams. Importantly, our analysis confirms the presence of finite size effects for small systems and sets a lower bound on the numbers of molecules that are needed to observe facsimiles of phase separation as diagnosed by the calculated pair distribution functions. The conclusions drawn from analysis of the pair distribution functions are reinforced in our analysis of the radial distribution function shown in S5 Fig.

Estimating the percolation transition line that delineates percolated and non-percolated networks

Associative polymers form networks characterized by physical crosslinks among stickers. Accordingly, we use the concept of a cluster, viz., a collection of unique chains connected via rotational interactions, to define the extent of percolation. In polymer melt simulations, the extent of percolation, known as the gel fraction in the polymer literature, is defined as the fraction of polymers participating in a percolating network that spans the simulation box in at least one direction [94]. More generally, we can use the fraction of polymers that make up the single largest cluster to quantify the onset of percolation and the changes to the extent of networking beyond the percolation threshold [95]. A molecular network cannot percolate the whole simulation cell when dilute and dense phases coexist. Accordingly, we choose the second definition for the order parameter that describes the percolation transition, and we denote this as ϕ_c [29].

Semenov and Rubinstein demonstrated that a percolation transition is purely a connectivity transition [45]. This implies that the identification of the percolation threshold is not achievable using a bona fide order parameter but instead relies on a suitable topological description. Here, we employ the midpoint of the ϕ_c vs. concentration curve to assess the onset of percolation and the percolation line or curve is obtained as the locus of points in the phase diagram for which ϕ_c = 0.5. In a system where finite size effects are minimized, the percolation transition is sharp having either a hyperbolic or sigmoidal shape as a function of concentration. Accordingly, the location of the percolation line will be relatively robust to the choice one makes for the percolation threshold.

A FUS-like system as an example of a linear multivalent protein

Wang et al. [24] recently uncovered the molecular grammar that contributes to the driving forces for phase separation of the protein FUS and a family of related proteins. They showed that, to first order, c_sat ∝(s_Ys_R)⁻¹, where c_sat is the measured saturation concentration of the FUS family proteins and s_Y and s_R are the number of tyrosine (Tyr) and arginine (Arg) residues, respectively. In FUS and other proteins of similar architectures, the Tyr residues are located primarily within the N-terminal disordered prion-like domain (PLD), whereas the Arg residues are located primarily within the partially disordered C-terminal RNA binding domain (RBD).

Using mutagenesis experiments, Wang et al. established that Tyr and Arg residues are the stickers in the FUS family proteins. Accordingly, the zeroth-order stickers and spacers representation used to model FUS in LASSI comprises of two parts: An N-terminal mimic of the PLD encompassing Tyr residues as stickers and a C-terminal mimic of the RBD that encompasses Arg residues as stickers. Wang et al. also measured c_sat for a 1:1 mixture of independent PLDs and RBDs interacting in trans. The c_sat for this system is approximately twice that of the c_sat for full-length FUS. Given the block copolymeric architecture of FUS, we denote the PLD and RBD as A_n and B_n, respectively for A and B-blocks of valence n. The model system of PLDs and RBDs interacting in trans is denoted as A_n+B_n (Fig 6A), whereas the system mimicking full-length FUS where the stickers can interact in cis and in trans is denoted as A_n-B_n (Fig 6B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Architecture of the linear multivalent systems.

(a, b) Cartoons to depict the A_n+B_n and A_n-B_n systems, respectively. Different colors of beads denote different species of stickers. Note that A_n-B_n can be simply considered as A_n+B_n where the two different sections of the proteins are joined together. (c) Linker lengths involved in the architecture (see also Table 1). Each sticker has a neighboring spacer bead that is 1 lattice site apart whereas the neighboring spacer beads are 4 lattice sites apart. This means that consecutive stickers are 6 lattice sites apart and also that the linkers connecting the two have a positive effective solvation volume.

https://doi.org/10.1371/journal.pcbi.1007028.g006

Within A_n and B_n blocks, spacers provide a uniform spacing of six lattice sites between stickers along the chain. We model spacers using a hybrid approach whereby a neutral spacer monomer is attached to each sticker with spacing of a single lattice site (Fig 6C). This choice was made to provide a distinction between A_n-B_n and A_n+B_n. Accordingly, the relative concentration of neutral beads will be higher in A_n-B_n when compared to A_n+B_n. This allows us to study linker-mediated differences between the driving forces for phase separation for A_n-B_n vs. A_n+B_n

Fig 7 shows phase diagrams for the A_n+B_n and A_n-B_n systems calculated using data from LASSI-based simulations. In panels (a) and (b), the ordinate quantifies the reduced temperature T* calculated as where ε is the effective energy of pairwise interactions between stickers from the A_n and B_n blocks. Panel (a) shows results for the A_n+B_n system. The bulk concentration in the A_n+B_n system is quantified along the abscissa as where are the bulk concentrations of A_n and B_n, respectively. Panel (b) shows the phase diagram for the A_n-B_n system where the abscissa represents the bulk concentration of this system.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Phase behavior of the linear multivalent systems.

(a, b) Phase diagrams for the A_n+B_n and A_n-B_n systems, respectively. The purple line is a 2-dimensional linear interpolation for = 0.025, and the area encapsulated by the purple line are where the systems have large density inhomogeneities and are thus considered to be phase separated. The green line is a 2-dimensional linear interpolation for ϕ_c = 0.5 and thus is the proxy for the percolation line. (c, d) and ϕ_c curves as a function of concentrations at T* = 0.383 (solid lines in (a) and (b)). (e) Width of the two-phase regime, w(T*), as a function of the reduced temperature. Not only does the A_n-B_n system have a higher critical temperature (T* ~ 0.6 vs. T* ~ 0.4), but also has a wider two-phase regime than the A_n+B_n system.

https://doi.org/10.1371/journal.pcbi.1007028.g007

Experiments show that the driving forces for phase separation are roughly twice as large for the full-length FUS compared to the system comprising of a 1:1 mixture of PLDs and RBDs [24]. This feature is recapitulated in LASSI simulations. For example, the width of the two-phase regime is larger for the A_n-B_n system compared to the A_n+B_n system for all values of T* as shown in panel (e) of Fig 7. The critical temperature is higher for the A_n-B_n vs. A_n+B_n system ( vs. , respectively). The valence of stickers is the main determinant of the concentration at the critical point whereas the interactions mediated by spacers determine the density inhomogeneities and the critical temperature. The impact of longer chains and increased valence of stickers per chain is also evident from the percolation threshold, which is crossed at a bulk protein concentration that is two-fold lower for the A_n-B_n system when compared to the A_n+B_n system, across all the simulation temperatures. Differences between the two systems are also evident in the degree of cooperativity of phase separation and the percolation transition as shown in panels (d) and (e) of Fig 7.

For each system, the intersection of the percolation threshold line with the two-phase regime shows that the dense phase predominantly forms a percolated droplet–panels (a) and (b) in Fig 7. Therefore, while phase separation without percolation is realizable, this is not the dominant scenario for associative polymers, where phase separation and percolation are typically conjoined to give rise to percolated droplets. The density of proteins in these percolated droplets is governed by the interaction strengths, modulated by T* and the effective solvation volumes of spacers [28, 29]. Unlike homopolymers, which comprise entirely of stickers or spacers depending on the solvent quality, associative polymers encompass a mixture of stickers and spacers. Stickers provide the driving forces for networking via reversible crosslinks and spacers determine whether or not these driving forces lead to phase separation via condensation. Indeed, the importance of sticker-driven percolation is evidenced in the persistence of percolated networks for both systems at high values of T*.

The observation that dense phases form percolated droplets has several implications: (1) on timescales that are concordant with or smaller than the average lifetime of physical crosslinks between stickers, the condensates will have elastic properties; this will be replaced by viscous behavior on timescales that are longer than the average lifetime of physical crosslinks [96]; (2) accordingly, condensates will have an intrinsic tendency for viscoelasticity [97] and long-lived crosslinks will cause hardening behavior as has been observed in many systems [1, 6, 7, 11, 22, 24, 39–41]; (3) the extent of crosslinking above the percolation threshold will change continuously with concentration [16, 29], and this will govern the overall structure, internal dynamics, and porosity of condensates; (4) reactions within condensates are likely to be constrained by the network topology formed as a result of inter-sticker interactions [98]; these constraints can create a variety of interesting kinetic signatures for reactions [99], including temporal memories as has been demonstrated recently for a system that undergoes thermoresponsive phase behavior [66]. Clearly, any description of biochemical reactions within condensates has to account for the structural and dynamical attributes of percolated droplets that are best described as network fluids.

Move set frequencies and diagnostics of converged simulations

We used results from simulations of linear multivalent protein system to assess the design of LASSI. The frequencies of the different move sets for simulations of the linear multivalent protein system are summarized in Table 2. Considerations that go into the design of move sets include the achievement of converged equilibrium distributions, with maximal computational efficiency, for each bulk concentration. Details of these considerations were described in the methods section. Diagnostics from short simulations are often useful to optimize the move set frequencies especially if multiple short trials are performed using very different starting configurations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Move frequencies according to their types.

They are normalized to the sum of all frequencies used in each simulation.

https://doi.org/10.1371/journal.pcbi.1007028.t002

Fig 8 shows the concentration dependence of acceptance ratios for each of the move set types, diagnosed for simulations of the A_n+B_n and A_n-B_n systems. The acceptance ratios show similar qualitative trends for both systems, even though there are clear quantitative differences. The move with the highest acceptance ratio in the dense regime is the double pivot move, signifying that the systems are transitioning into a pure polymer melt. The second most accepted move is the local move; extrapolating from the higher concentrations it is expected that the acceptance of local moves should also become small and that the double pivot move will be the most dominant way to alter chain configurations, since even the move of an individual monomer will require that multiple monomers from multiple chains are moved simultaneously. Both systems have similar qualitative trends for the translation move where we see a transition from being accepted at low concentrations to not being accepted at higher concentrations. Since the proteins in the A_n-B_n system are twice long as the A_n+B_n system, the absolute acceptance ratio of the translation move is always lower in the A_n-B_n system.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Analysis of acceptance ratios for different move sets.

Curves with different colors indicate acceptance ratios of different types of moves. The dashed lines show the saturation concentrations. The data are obtained from simulations with T* = 0.383. (a) Acceptance ratio data for the A_n+B_n system. (b) Acceptance ratio data for the A_n-B_n system.

https://doi.org/10.1371/journal.pcbi.1007028.g008

Analysis of acceptance ratios of different move sets within droplets will be helpful for estimating correlation lengths and amplitudes of conformational and concentration fluctuations within droplets. Cluster moves have high acceptance ratios in the dilute regime whereas the acceptance ratio nearly vanishes as the concentration increases. This is intuitive since the likelihood of steric clashes increases with a decrease in available volume and this is coupled to the simultaneous increase in the fraction of molecules in the largest cluster. We note here that the cluster moves have the most dramatic change in acceptance ratios from values near 1 to values near 0. However, the apparent inefficiency of cluster moves in dense configurations cannot be used as a rationale to quench such moves. In fact, as shown in panel (a) of Fig 9, phase separation, diagnosed in terms of , is suppressed if cluster moves are not part of the move set. This highlights the importance of cluster moves for generating bona fide phase separation as these facilitate coalescence that leads to condensation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Importance of cluster translation moves.

(a) and (b) ϕ_c curves for the A_n+B_n (purple) and A_n-B_n systems (green) at T* = 0.383. The solid lines are identical with the curves in panels (c) and (d) of Fig 7. The dotted lines show the simulation results under the same system conditions but the frequency for cluster translation moves is set to zero. Not only do the systems phase separate and percolate at higher saturation concentrations, but also we can see that both percolation and separation are suppressed highly. Furthermore, errors are generally higher, due to the systems being highly dependent on the initial conditions of the system.

https://doi.org/10.1371/journal.pcbi.1007028.g009

Mixtures of N130 and the rpL5 peptide as an example of a branched multivalent protein system

LASSI can also be deployed to study the phase behavior of branched multivalent proteins that undergo phase separation via obligate heterotypic interactions. Examples of branched multivalent proteins are molecules that form symmetric, stable oligomers such as nucleophosmin 1 (NPM1) and synthetic systems such as the corelets designed by Bracha et al. [100]. NPM1 is a key molecule within the granular component (GC) of the nucleolus [101]. Three coexisting layers define the nucleolus and the GC is the outermost layer. Within the GC, NPM1 appears to form facsimiles of percolated droplets in complex ribosomal proteins with Arg-rich motifs [17, 30]. A minimalist system that mimics the phase behavior of the GC comprises of a truncated version of NPM1, referred to as N130, and an Arg-rich peptide, designated as rpL5 [31–33]. Both NPM1 and N130 form symmetric pentamers in the presence of Arg-rich peptides [102]. The pentamer formed by the association of folded domains serves as a scaffold for displaying disordered C-terminal tails that are defined by two distinct acidic tracts. The system also features an N-terminal disordered region with a well-defined acidic motif.

The FUS system is an example of a protein that undergoes phase separation via obligate homotypic interactions. This implies that the interactions necessary and sufficient for driving phase separation are encoded within the sequence of FUS and the strengths of these interactions can be modulated by changes to solution conditions. The N130 + rpL5 system is an example of a system that undergoes phase separation mainly via obligate heterotypic interactions that involve interactions between residues in the acidic tracts of N130 and the Arg motifs of rpL5. This could be viewed as an example of phase separation via complex coacervation, providing the heterotypic interactions are purely electrostatic in nature [31–33]. However, in general, there is likely to be combination of long- and short-range interactions that contribute to the spectrum of heterotypic interactions, and hence we refer to this class of molecules as drivers of phase separation via obligate heterotypic interactions.

In addition to demonstrating the applicability of the LASSI engine for simulations of branched systems, we use the analysis as an opportunity to highlight key conceptual features of multicomponent systems that undergo phase separation via obligate heterotypic interactions. There are three main features that are borne out in the analysis: (1) For fixed temperature, the order parameters are the concentrations of the proteins that drive phase separation via obligate heterotypic interactions. In such systems, the coexistence curves, i.e., the binodals, will have a closed loop form. These will be ellipses for two components and n-ellipsoids for systems that involve up to n obligate heterotypic interactions to drive phase separation. (2) The systems will support reentrant phase behavior as has been reported recently for protein-RNA mixtures that undergo phase separation via obligate heterotypic interactions [103]. (3) The apparent saturation concentration of a component molecule in a system that undergoes phase separation via obligate heterotypic interactions will show non-trivial dependencies on its bulk concentration. These dependencies are governed directly by the slopes of the tie lines that pass through the point corresponding to the bulk concentration and intersect the binodal at coexisting concentrations that equalize the chemical potentials of all species in the dense and dilute phases. Here, we use the example of the N130 + rpL5 system to showcase the three central features of phase diagrams for systems that undergo phase separation via obligate heterotypic interactions.

In the LASSI representation, N130 pentamers with disordered tails are modeled using a five-armed structure. This approach follows the strategy of Feric et al. [30], which was based on the fact that pentamers do not dissociate under conditions where NPM1 / N130 undergo phase separation. Each arm comprises of two sticker sites to mimic the presence of the A1 and A2 acidic tracts within the disordered tails of NPM1 / N130. Therefore, each N130 pentamer displays a total of ten sticker sites. The spacers between each A1 tract and the N130 core as well as between each pair of A1 and A2 tracts on a disordered tail are phantom tethers, which means that their effective solvation volumes [29] are set to zero. Each rpL5 peptide has two sticker sites corresponding to the two Arg-rich motifs along the sequence. Schematic representations of the coarse-grained architecture used for N130 and rpL5 are shown in Fig 10, and the move set frequencies are summarized in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Architecture of an archetypal branched multivalent system.

(a) Schematic to depict the overall architecture. The pentamer with 10 orange stickers represents the N130 molecule where the gray central oligomerization domain is modeled as a neutral spacer monomer, and the rpL5 peptide is modeled as a linear molecule with 2 blue stickers. (b) Relevant length scales for the architecture (see also Table 1). For the rpL5 molecule a linker length of 3 was chosen between the two stickers, and for the N130 molecule the first sticker (modeling the A1 tract) is 1 lattice site away from the hub spacer whereas the second sticker (modeling the A2 tract) is 3 lattice sites away from the first sticker.

https://doi.org/10.1371/journal.pcbi.1007028.g010

Percolation lines have parabolic shapes

The percolation line, constructed as a function of the concentrations of two multivalent molecules, has an overall parabolic shape (panel (a) of Fig 11). This feature may be explained as follows: Let f₁ and f₂ denote the fractions of N130 and rpL5 molecules that are bound in a network; s₁ and s₂ will denote the valence of stickers on N130 and rpL5, respectively; for the current implementation of the N130 + rpL5 system, s₁ = 10 and s₂ = 2. Based on the Flory-Stockmayer theory [49, 51], the percolation threshold is crossed when f₁f₂(s₁−1)(s₂−1) > 1. If we keep (s₁−1)(s₂−1) constant, the threshold concentration for percolation will be governed by the product f₁f₂. Accordingly, if there is a large excess of N130 (component 1) compared to rpL5 (component 2), then f₁ → 0 and f₂ → 1, and the system does not undergo a percolation transition. In this scenario, every rpL5 molecule is crosslinked to two sticker sites from one or two N130 molecules. However, since the relative stoichiometry favors N130 molecules, there is a vast excess of un-crosslinked N130 molecules and the network cannot grow. Percolation is also inhibited when the converse situation arises, i.e., when there is a large excess of component 2 with respect to component 1. Accordingly, the percolation line takes on a parabolic form in the plane defined by the concentrations of N130 and rpL5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Phase behaviors of the branched multivalent systems for T* = 0.25.

(a) Full phase diagram, where the purple line denotes the proxy for the binodal and the green line is the proxy for the percolation line (see also the caption for Fig 7). The phase-separated region has an elliptical shape and we have a closed loop, which demonstrates re-entrant phase behavior, whereas the percolation line has a conical shape extending into much higher densities. The solid black lines denote contours of constant total concentration where L1 is the lowest concentration and L3 is the highest concentration. Note that both axes are represented in the log scale. (b, c) and ϕ_c curves as a function of relative stoichiometric ratio of N130 and rpL5 along the constant-concentration contours. (d) Plot of Λ vs. the apparent stoichiometry along lines L1, L2, and L3.

https://doi.org/10.1371/journal.pcbi.1007028.g011

Binodals for systems defined by obligate heterotypic interactions will form closed loop ellipses

Given that the phase behavior of the N130 + rpL5 system is driven by heterotypic interactions involving the A1 / A2 tracts from the N130 tails and the Arg-motifs from rpL5, we constructed binodals by keeping the simulation temperature fixed and varied the concentrations of N130 and rpL5 molecules. Phase diagrams defined by N130 concentration along the abscissa and rpL5 concentration along the ordinate are shown in panel (a) of Fig 11. The general shape of the binodal is comparable with that of the experimentally determined phase diagram [33], even though direct comparison is not straightforward because the scarcity of experimental data points does not yield a full binodal.

The phase boundary, defined by the density transition, is an ellipse that forms a closed loop in the plane defined by the concentrations c₁ and c₂ of N130 and rpL5, respectively. In associative polymers, the phase behavior is governed by the affinity between stickers, the valence of stickers, and the effective solvation volumes of spacers [28, 29]. For fixed c₁ that intersects the two-phase regime an increase in c₂ will lead to an entry into the two-phase regime caused by a density transition as c₂ approaches c₁. However, as c₂ increases well beyond c₁, the joint system exits the two-phase regime. This is because phase separation is driven by obligate heterotypic interactions and while there is a growing excess of rpL5 molecules there are not enough N130 molecules to drive the density transition via inter-sticker interactions. Similar reasoning applies to describe the reentrant behavior that will result by keeping c₂ fixed at a value that intersects the two-phase regime and increasing c₁.

Taken together, the parabolic percolation lines and elliptic forms for two-phase regimes define conic sections that highlight reentrant phase behavior whereby fixing the concentration of component 1 and increasing the concentration of the second species can lead to phase separation and percolation at a low threshold concentration of component 2 and exit into the one-phase, non-percolated regime beyond a second higher threshold concentration for component 2. This type of reentrant phase behavior, will be a general feature of multicomponent systems that undergo phase separation via obligate heterotypic interactions; indeed, reentrant phase behavior has been reported for a model protein + RNA system [103].

Apparent stoichiometric ratios can be different from effective stoichiometric ratios

Stoichiometry of molecules that drive phase separation is another key parameter that determines the functions of biomolecular condensates formed by multicomponent systems [104]. The apparent stoichiometric ratio is calculated as the ratio of the concentrations of stickers of types s₁ and s₂ for N130 and rpL5, respectively such that . However, the effective stoichiometric ratio can be different from if excluded volume effects modulate the effective concentration of stickers. We fit an ellipse to the two-phase boundary and determined the major axis of this ellipse. The effective stoichiometric ratio should be unity along the major axis. As shown in panel (a) of Fig 11, the major axis deviates from the line along which . Therefore, and angle between the major axis and the line along which quantifies the impact of excluded volume on changes to effective concentrations of stickers that in turn modifies the stoichiometric ratios.

The synergy between stoichiometry and phase behavior can be analyzed by quantifying the order parameter and the topological parameter ϕ_c as a function of apparent stoichiometry for fixed bulk concentration. Along each gray line in panel (a) of Fig 11 the total concentration defined as c_bulk = (c₁c₂)^½ is fixed, although the stoichiometries will vary. The mean values of c_bulk along L1, L2, and L3 are 2.09×10⁻² (voxel^–1), 2.46×10⁻³ (voxel^–1), and 3.33×10⁻⁴ (voxel^–1), respectively and the value of ranges from 0.36 to 22.62 along each of L1, L2, and L3. Panels (b) and (c) in Fig 11 show the variation of and ϕ_c as increases along L1, L2, and L3, respectively. Along L1, the value of is essentially zero irrespective of stoichiometry because L1 lies is outside the two-phase regime. However, a system spanning percolated network forms for stoichiometries in the range 1.2 ≤ ν₁₂ ≤ 13 along L1. This is because the concentrations of both components are well above the percolation threshold along L1 thus ensuring that stickers readily find one another even without a density transition. In direct contrast, along L3, we observe phase separation, characterized by values of > 0.025 for a range of stoichiometries, but none of these support the formation of a percolated droplet (ϕ_c < 0.5 for all stoichiometries). Along L2, we observe phase separation for stoichiometries in the range 1.15 ≤ ν₁₂ ≤ 16 and percolation for stoichiometries in the range 2.14 ≤ ν₁₂ ≤ 11.3 such that phase separation enables the formation of a percolated droplet.

In panel (d) of Fig 11, we introduce a new structural parameter Λ, which we define as a convolution of and ϕ_c such that . Here, the convolution is calculated as a logical AND gate, which becomes a simple product. The parameter Λ quantifies the convolution of the density and network transition and provides an estimate of the extent to which the phase separation and percolation are coupled as the apparent stoichiometry is varied for a fixed bulk concentration. The profile of Λ is reminiscent of profiles measured by Case et al. [104] for the dwell time of signaling molecules as a function of stoichiometric ratios that govern the formation of condensates at membranes. This suggests that dwell times, which are experimentally accessible parameters, might actually be proxies for the structural features of the condensates as measured by the convolution between phase separation and percolation and the extent of network formation within the condensate.

The key finding is that the combination of the bulk concentration and stoichiometric ratio (as opposed to stoichiometry alone) will determine the quench depth into the two-phase regime. This in turn determines whether a system-spanning network forms without phase separation or if phase separation enables the formation of a droplet-spanning network. The structure of condensates and the overall phase behavior cannot be fully described in terms of c_bulk or ν₁₂ alone. Instead, this requires the consideration of both parameters jointly and relative to the quench depth, which refers to the location in the two-phase regime and with respect to the percolation line. This is important because the extent of crosslinking and the time scales associated with crosslinks will determine the material properties of the condensate. This in turn should contribute to parameters such as the dwell times of clients within condensates [104].

Saturation concentrations need not be fixed parameters in multicomponent systems

The concept of a saturation concentration is one of the defining hallmarks of phase separation [2, 27]. For fixed solution conditions, phase separation in a closed two-component system (or pseudo one-component system) comprising of a protein and solvent is realized when the bulk concentration of the protein denoted as c exceeds a saturation concentration denoted as c_sat. The presence of a saturation concentration can be quantified by measuring the concentration of protein in the coexisting dilute phase as c increases. This value will not go above c_sat. Strikingly, the presence of a saturation concentration has been confirmed in living cells for model disordered proteins by Brangwynne and coworkers using optogenetic tools in mammalian cells [100, 105] and by Khan et al. [106] using yeast as a model system.

If the protein whose phase behavior is being interrogated is part of a system where obligate heterotypic interactions drive phase separation, then whether or not the concept of a saturation concentration continues to be valid will depend on the nature of the binodals. We illustrate this by coopting the elliptic phase boundary from Fig 11 for a three-component system that comprises of N130, (component A), rpL5 (component B), plus a solvent that is implicit in the LASSI simulations. This 3-component system may be thought of as a pseudo two-component system. For fixed temperature, the top row of panels (a)-(c) in Fig 12 show three types of elliptical, closed loop binodals. These are constructed in a plane where [B] increases along the positive direction of the abscissa and [A] increases along the positive direction of the ordinate. The bottom row in each panel shows the result to be expected were we to measure the concentration of A in the dilute phase, designated as [A]_Sol, as the bulk concentration [A] is varied. In each of these measurements, the concentration of B is fixed at a specific value.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Slopes of tie lines within elliptic binodals are important for systems that undergo phase separation via obligate heterotypic interactions.

The ellipse is drawn to fit the locus of points based on the value that meets our criteria for a density transition (see main text). Data for constructing the ellipse were taken from simulations of the N130 + rpL5 system–see Fig 11(A). This ellipse is used to assess the impact of slopes of tie lines for a two-component system comprising of macromolecule A that undergoes phase separation via obligate heterotypic interactions with macromolecule B. (a) Ellipse with nearly horizontal tie lines. The vertical lines shown in green, grey, and blue correspond to fixed values for [B] along the abscissa. As [A] increases, the system traverses across the two-phase regime, delineated by the ellipse, starting outside the ellipse, crossing the ellipse, and exiting the ellipse at high concentrations of A. (b) For each fixed value of [A], the plot shows how [A]_Sol varies with [A]. The red points on each curve were extracted from within the two-phase regime, whereas the black points lie outside the two-phase regime. Clearly, [A]_Sol does not stay fixed as [A] increases. (c) Equivalent plot to that shown in panel (a) for the tie lines that we obtain for the N130 + rpL5 system. (d) Equivalent plot to that shown in panel (b). Note the non-linear variation of [A]_Sol as [A] increases. (e) Ellipse annotated with vertical tie lines. In this case, phase separation of [A] does not depend on obligate heterotypic interactions with B, but B can bind to A and has a choice of binding preferentially to A in either the dense or dilute phase. Here, A becomes the macromolecule and B the ligand. (f) Preferential binding of the ligand to the macromolecule in its dilute phase will shift the saturation concentration, assessed in terms of [A]_Sol, upward and this shift will depend on [B]. Accordingly, the plateau value of [A]_Sol in the two-phase regime shifts to higher values for higher values of [B].

https://doi.org/10.1371/journal.pcbi.1007028.g012

Slopes of tie lines within the elliptic binodals determine how [A]_Sol varies with [A] in the two-phase regime

For concentrations of A and B that place the pseudo two-component system in the two-phase regime–red points along each of the curves in the bottom rows of panels (a)–(c)–we find that [A]_Sol can change as [A] increases. If the tie lines are horizontal or nearly horizontal, then [A]_Sol will vary linearly with [A]. Non-linear variations of [A]_Sol with [A] will result for tie lines with positive or negative slopes. This is shown in panel (b) for tie lines with positive slopes. If the tie lines are essentially vertical, then the standard expectation regarding the invariance of [A]_Sol with [A] within the two-phase regime is recovered. However, even in this scenario, the plateau value of [A]_Sol will shift upward or downward as the value of [B] increases–the upward shift is shown in panel (c) of Fig 12. Here, B acts as a bona fide ligand for A, which is the macromolecule. Preferential binding of B to the dilute phase leads to an increase in [A]_Sol as depicted in panel (c) of Fig 12. Ligand-mediated shifts in saturation concentrations arise due to polyphasic linkage, a phenomenon first introduced by Wyman and Gill [107].

The main conclusion is that the concept of a saturation concentration, as defined for a pseudo one-component system, does not transfer over to multicomponent systems where phase transitions are driven by obligate heterotypic interactions. Instead, the slopes of tie lines or the geometries of tie planes in higher dimensional ellipsoids will have a direct bearing on inferences from measurements where the bulk concentration of a protein or RNA component is varied when condensates are observed and the concentration of the molecule of interest is quantified in the coexisting dilute phase. This insight emerges from our ability to deploy LASSI to compute full binodals for multicomponent systems.