■ INTRODUCTION

Intrinsically disordered proteins (IDP) are able to spontaneously self-assemble into membraneless organelles or droplets which exhibit liquid-like behavior. Several studies indicate that these membraneless assemblies are formed by a liquid-liquid phase separation (LLPS) process. [1][2][3][4][5][6][7] Liquid droplets formed by IDPs dynamically exchange molecules with the surrounding environment, amalgamate into larger droplets, and can be deformed by flows. 4,6,8,9 These droplets are known to become more viscoelastic, transform into a gel-like state with time, and eventually behave as solids. 4,10 The current consensus in the literature is that phase separation in biomolecular systems is unlike the conventional LLPS observed in solutions of uncharged polymers. 11 The strong interactions among multiple components in the biomolecular systems lead to more complex phase diagrams. 7,[12][13][14][15] Analogous to the well-known phase behavior in solutions of associating polymers, [16][17][18][19] the role of percolation in concentrated solutions of biomolecules has been recognized to be an important factor to understand the phase behavior of biomolecular systems. Modeling of biomolecules as simple sticker-spacer chains has revealed two possible mechanisms for the formation of biomolecular condensates: system-spanning percolated network formation without any phase separation, and phase separation driven percolation with multivalent proteins condensing into dense droplets with internal percolated structure. 14,[20][21][22] Numerical studies have also uncovered the possibility of metastable microcluster phases in these multicomponent protein systems at low to intermediate concentrations. 23 In general, the phase behavior of IDP solutions arises from a confluence of their several key attributes. These include electrostatic and hydrophobic interactions among the residues, chemical sequence, and the extent of blockiness of sequences of a particular kind (positively charged, negatively charged, hydrogen-bonded, and hydrophobic). As a result of this confluence, the phase behavior can be quite complex in realistic experimental situations. Even for the simplest situation of solutions of uncharged flexible homopolymers, such as polystyrene in carbon disulfide, 16 the phase diagram is as sketched in Figure 1(a) as a plot of temperature versus polymer concentration. In this simple version of phase behavior of associating polymers, the phase diagram displays a homogeneous solution phase, percolated gel phase, percolation line separating solution and gel phases, and the coexistence of a dilute solution phase with a gel phase at lower temperatures. The meeting point of the percolation line at the critical point for coexistence of the solution and gel phases is a tricritical point. Inside the coexistence curve, the equilibrium situation is the coexistence of one solution phase and one gel phase at a fixed temperature. In reaching this equilibrium, it is not uncommon that the kinetic process of phase separation into the solution and gel phases can be slow due to barriers arising from the necessity of structural reorganization to reach equilibrium.

In the presence of charged residues, association of oppositely charged domains either from the same molecule or from several interpenetrating molecules can readily occur, as sketched in Figure 1(b), by releasing their relevant counterions. 24 Such associations formed by multiple ion-pairs can spontaneously result in formation of networks involving many chains. 24 Depending on the availability of such molecules in the system, these associations can lead to either clusters of finite size or a percolating gel. Furthermore, if the sequence of the molecule contains domains of different kinds, then a collection of large number of such molecules can exhibit the well-known phenomenon of microphase separation. 26,27 Even in the simplest situation of diblock copolymers (one block is charged with a given sign, and the other block is uncharged), the system exhibits well-organized structures, without the LLPS (macrophase separation). 25 The canonical morphologies formed by diblock copolymers with little solvent are lamellae, cylinders, and spheres, as sketched in Figure 1(c). The signature of microphase separation is that the scattering intensity shows a peak at a finite scattering wave vector q corresponding to a distance scale comparable to the radius of gyration of the molecule. This is in contrast with the LLPS (macrophase separation), where the scattering intensity does not show any peak at finite q and it decays monotonically with q from its value at q = 0, representing that the resultant structures are of macroscopic length scale. In the presence of increasing amount of solvent, the first order transition line between the disordered and microphase separated microstructures reaches the macrophase phase separated critical point at the Lifshitz point. Thus, in general, the phase behavior of IDP solutions can be very complex because all of the issues depicted in Figure 1 are simultaneously operative.

The primary goal of the present paper is to address the contributions of the various concepts underlying the sketches in Figure 1 to the formation of emergent structure from IDP molecules. Toward this goal, we model finite-sized droplets of varying IDP concentration and account for hydrophobic and electrostatic interactions, and sequence, to obtain the eventual structural organization of the molecules inside the droplets. The task of determination of the full phase diagram is more challenging and is relegated to future work. Already with the goal of discerning structure of interpenetrating IDP, even in the absence of other cofactors or crowding agents, we find emergence of microphase separated microstructures. This new feature and the various ingredients mentioned in Figure 1 need to be accounted for in advancing an understanding of phase behavior of IDP solutions.

In our investigation, we have chosen α-synuclein as the IDP, because it is a natively unstructured protein 28,29 and is the prime example of IDPs. Furthermore, α-synuclein is one of the well investigated systems involving IDPs. A brief summary of the literature on α-synuclein is as follows. It is known to misfold and aggregate into insoluble, highly structured, and βsheet-containing cytotoxic oligomers and amyloid fibrils which form Lewy bodies and Lewy neurites, the hallmarks of Parkinson's disease (PD) and other synucleinopathies. [30][31][32] In monomeric form, it consists of three primary domains; Nterminal region with overall positive charge, the hydrophobic nonamyloid-β component (NAC) region and the C-terminus, a mostly negatively charged region. Due to highly concentrated charges at the termini and high overall hydrophobicity of the NAC region, α-synuclein has the characteristics of an IDP. The transient and dynamic electrostatic and hydrophobic intramolecular interactions are significantly influenced by the surrounding environment and lead to an ensemble of monomeric conformational states. This in turn influences the aggregation pathways and whether these are potentially neurotoxic or neuroprotective. 9,31,[33][34][35][36] Determining the conformations within the aggregates or the conditions and cofactors that can destabilize these aggregates will aid in the study of antiaggregation therapeutics.

A large body of studies suggest that α-synuclein undergoes LLPS by forming liquid droplets similar to many other IDPs at the early stages, which mature into hydrogels containing fibrillar aggregates. 4,9,10,36 Recent experimental studies do suggest that at physiological pH (∼7.4) α-synuclein is highly soluble and can undergo phase separation only in the presence of additional factors decreasing the free energy barrier for LLPS. 4,9 Consequently, droplets in bare α-synuclein solutions without any crowding agents are observed only at relatively higher concentrations. 4 Further, fibrillation of α-synuclein is also pH dependent and occurs at increasing rates at acidic pH as is known from in vitro experiments. 37,38 The acidic Cterminus residues are largely negatively charged at neutral pH and become charge neutral at low pH. Such a reduction in charge density reduces intra-and intermolecular electrostatic interactions and enhances hydrophobicity responsible for fibrillation. α-Synuclein's involvement in the pathogenesis of PD has led to significant experimental and theoretical efforts to understand the aggregation mechanism of α-synuclein. Yet, complete understanding of the mechanisms that result in aggregation into droplets, and their phase behavior remains still unclear and needs further extensive research.

In the present study, we consider α-synuclein solution at physiological conditions and investigate the self-assembly behavior inside a finite volume using extensive multiscale numerical simulations. We observe that at these conditions, electrostatically driven transient microclusters emerge if the protein concentration is low, and percolated structures consisting of microdomains emerge if the protein concentration is high. Moreover, varying salt concentration reveals that these microstructures become looser with increasing ionic strength of the solution. With our primary focus being the structural elucidation of the protein-dense phase, we show below that the dense phase exhibits microphase separation. Hence the feature of microphase separation needs to be integrated with studies of phase diagrams of IDPs and other biomacromolecular systems.

■ NUMERICAL MODEL

We model the protein as an united atom polymer, where each amino acid residue is represented by one single bead and takes the position of C α . For simplicity, we consider each bead with the same diameter σ = 0.38 nm and same mass m = 103.29 g/ mol such that the molecular weight of the chain is equal to that of α-synuclein (14460 g/mol). Our united atom model incorporates force fields to represent bonding, backbone rigidity, and nonbonded hydrophobic and electrostatic interactions. The bonded interactions are implemented by the harmonic potential

where r is the distance between two neighboring beads along the chain. The bond constant K is taken to be 121.11 kcal nm -2 mol -1 and the equilibrium bond length r 0 is equal to the bead diameter σ. The bending energy and the torsional energy of a chain is modeled by sequence-specific statistical potentials derived from the coil library. 39 Ramachandran plots are extracted from the random coil database and we use them as inputs to generate the coarse grained bending and torsional potentials as described in ref 40. In particular, we adopt triple combinations of Glycine (G), Proline (P), and "Generic" (X) (referring to the remaining 18 amino acids) for neighboring 3residue fragment forming a bending angle to describe the corresponding potential. The torsion potentials are considered for double combinations of amino-acid residues in the middle of a 4-residue fragment O-Z-Z′-O forming the dihedral angle with {Z,Z′} ∈ {X,G,P} and {O,O′} are any of the amino acid residues. Electrostatic interactions are modeled using the Debye-Huckel potential, which takes into account the electrostatic screening due to the presence of salt, having the functional form (in units of k B T)

Here, ξ D is the Debye screening length and B = 0.7 nm is the Bjerrum length, corresponding to the dielectric constant of 80 for the solvent medium (water) at room temperature. q i and q j are the charges of beads i and j. We consider the standard charges at physiological pH with aspartic acid and glutamic acid having a charge q = -1, lysine and arginine with a charge q = +1, and histidine with q = +0.5. The electrostatic interaction is cutoff at a distance which varies depending on ξ D . In particular, we use a cutoff distance of 4.1 nm (ξ D = 1.36 nm), 3.61 nm (ξ D = 1 nm), 3.12 nm (ξ D = 0.785 nm), and 2.66 nm (ξ D = 0.68 nm) for the ionic strengths 50, 100, 150, and 200 mM, respectively. We model the nonbonded hydrophobic interactions using the potential

where U LJ is the standard Lennard-Jones potential given by

The energy parameter ε (in units of k B T) and dimensionless λ are chosen according to the Kim and Hummer (KH) model parametrized to describe protein-protein interactions by fitting to experimental data for osmotic second virial coefficient of lysozyme and the binding affinity of the ubiquitin-CUE complex. 41 In particular, we use the parameter set D of KH model as experiments suggest high solubility of α-synuclein at neutral pH and it has also been used earlier to appropriately describe the interactions between IDPs where all residues are essentially fully exposed to the solvent. 2 The cutoff distance for the hydrophobic interaction is 2.5σ beyond which it becomes negligible.

Each monomer bead of a protein chain undergoes Langevin dynamics at a temperature T according to the equation of motion

(5) where, γ denotes the friction coefficient, and -∇U indicates the forces arising from all the potentials introduced above. Following fluctuation-dissipation theorem, thermal fluctuations are incorporated by the Gaussian white-noise vector Γ, with ⟨Γ⟩ = 0 and

. We integrate the Langevin equation using a stochastic Verlettype algorithm at a constant temperature T = 300 K with a 20 fs time step and 0.4 ps -1 friction coefficient. 42 We characterized and validated the conformational properties of an isolated α-synuclein chain using our model (see Supporting Information, SI). Our choice of force-fields and parameters leads to good agreement with single-chain experiments in literature and are able to accurately capture the intrachain interactions. In all of our bulk simulations, while considering different protein concentrations, we keep the number of protein chains fixed at N = 1925 and vary the cubic simulation box size from 68.4 to 136.8 nm. For each parameter set considered to obtain the characteristic results, we take into account time average of the data and also ensemble average over 8-10 independent runs.

■ RESULTS AND DISCUSSION

Behavior. We have considered α-synuclein solutions with five different protein concentrations in the range C p ∈ (1.25, 10) mM. We can identify three scenarios in the protein solutions, a dispersed fluid-like structure at lower concentration, transient micro clusters at intermediate concentrations, and a large scale percolated structure at the highest concentration. Figure 2 shows the typical equilibrium snapshots of our simulations for the protein concentrations 1.25, 4.49, and 10 mM, where the protein chains have a color according to the cluster they belong to. We consider two protein chains to be part of the same cluster if any of the monomer beads from those two chains are within a cutoff distance 0.464 nm, which is the distance corresponding to the first minimum in the interchain radial distribution function. By performing cluster analysis with this criterion, each cluster can be characterized in terms of the cluster size (s), which is the number of protein chains belonging to that cluster. Each cluster in Figure 2 is colored in a RGB color-scale according to the sequential cluster index C id (s) after the clusters are sorted in a decreasing order of cluster size. The cluster consisting of maximum number of protein chains (s max ) with an index C id (s max ) = 1 is colored red and any cluster with s < 11 is colored blue.

We have obtained a quantitative understanding of the role of protein concentration on the aggregation-percolation phenomenon (Figure 2) from the cluster size distribution, P(s), where s is normalized by the total number of chains in the system and is presented in terms of the percentage of polymer chains belonging to a particular cluster. Figure 3 shows the cluster size distribution for five different protein concentrations at a fixed salt concentration, C s = 100 mM (∼ physiological concentration). It is evident that at C p = 10 mM in addition to small clusters a distinct broad peak appears at very large cluster size which is not present at lower concentrations. It reveals that almost 92-97% of protein chains are interconnected and belong to the same percolated large cluster. At the same time rest of the protein chains remain mostly in single chain form as can be observed from the inset plot. The inset plot shows the distribution as a function of the absolute number of chains belonging to a cluster instead of the percentage of chains and is zoomed into the lower cluster size region of the main plot. As the concentration increases from C p = 1.25 mM to 5.24 mM (lowest to intermediate), the distribution curves extends to higher cluster sizes with a slower decay of the curves and a slight nonmonotonicity at larger s-values becomes apparent, which is indicative of the microclusters at intermediate concentrations. Overall, we observe the tendency of microclusters emerging with increasing protein concentration, which transitions into a system-spanning percolated phase beyond a threshold concentration.

In order to get an insight into the aggregation mechanism, we analyzed the details of the structural organization of the protein chains at different concentrations of the α-synuclein solution. We characterize the structural correlations in our system in terms of the static structure factor, S(q) and the radial distribution function, g(r). For a set of N residues, the structure factor is defined as

where q is the scattering wave vector and r j , r k are the position vectors of the residues j and k, respectively. In our analysis, the interchain radial distribution function (RDF) is defined as

i j i j inter (7) where ρ(r) is the density inside a spherical shell with radius r and width 0.02 nm for an isotropic system. α β r i j is the distance between residues i α in chain α and j β in chain β. In Figure 4(a), we first determine the interchain radial distribution function of the protein chains to resolve the effect of concentration variation on structural correlation. We can observe a characteristic broad peak around 1 nm with decreasing height as the concentration increases. Interestingly, the fluid-like phase at the lowest concentration exhibits a stronger shortrange arrangement than the microclusters at intermediate concentrations and the percolated network at the highest concentration. It is evident from the cluster size distribution that the system consists of more dimeric and trimeric protein aggregates at lower concentrations. Protein chains in these aggregates are more compact and closely bound to each other as is evident from Figures S5,S6, and S7. With increasing protein concentration, hairpin-like configurations of proteins with intrachain contacts between the N and C-terminals are less prevalent. These structures acquire a slightly more extended configuration due to interpenetration by other chains and an increasing interchain repulsion from the surrounding Cterminals. This manifests in terms of a broad prominent peak with increased height for a protein residue after the small first peak as the concentration of proteins is lowered. In the fluidlike phase with a simple exponential fit starting at r = 2 nm, the structural correlation length is 2.64 nm for C p = 1.25 mM and 3.3 nm for C p = 2.96 mM. Figure 4(b) shows the structure factor of the system with protein concentration of C p = 10 mM to recognize the characteristic length scales of structural correlation in the percolated phase (see Figure S8 for other concentrations). We observe a prominent peak at q = 0.175 nm -1 and some additional small peaks toward higher q values. The first peak corresponds to an average distance spacing of 6 nm, which is almost 2 times the radius of gyration value R g = 2.85 nm of one single protein chain (see Table S1). Presence of the peaks at nonzero q-values indicates microstructural organization of the protein chains within the percolated phase.

Given the above observations of structural organization, we considered the contact map of all the 140 residues of αsynuclein as shown in Figure 5(a) for the system with C p = 10 mM and C s = 100 mM to characterize the relative arrangement of different residues, as it represents the likelihood of any two residues being close to each other. Here, two residues are regarded to be in contact if they are within a cutoff distance of 0.95 nm (the cutoff considered for the nonbonded hydrophobic interactions (U nb ) and roughly the second minimum position in interchain radial distribution function) and we have not accounted the contacts from the residues belonging to the same protein chain, i.e., only interchain contacts are considered. The color scale is according to the contact fraction, f c , which is the ratio of a particular residue-residue contact to that of the maximum residual contact (the highest residue-residue contact number out of all the combinations of residue pair). It is evident that the large scale percolated structure is formed mainly because of the contacts between the charged N-terminal and C-terminal regions of the protein chains. In addition, there are non-negligible number of contacts between the residues belonging to the hydrophobic NAC-region. The interchain contacts become more clear from the partial radial distribution function plot in Figure 5(b), where we consider interchain pair correlation between the three N-terminal, NAC, and C-terminal segments. As expected, we can observe the N and C-terminal radial distribution curve having the maximum peak around the distance 1 nm. This peak height decreases and also the 3-3 curve goes to the limit 1 faster with an increase in protein concentration as we show in Figure S7. This is in accordance with Figure 4(a) and our discussion that repulsion from C-terminal results in structural rearrangement at higher concentrations.

Effects of Salt. So far, we have considered a fixed salt concentration of C s = 100 mM for the systems with different protein concentrations. To elucidate the effect of salt, we study the structural properties under the variation of salt concentration in the solution while keeping the protein concentration fixed at C p = 10 mM. The cluster size distributions of the protein chains for the salt concentrations C s = 50, 100, 150, and 200 mM are shown in Figure 6(a). Here, again, we normalize the cluster size s by the total number of chains in the system and is presented in terms of chain percentage. There is a distinct peak at higher chain fraction for clusters containing more than 90% of the protein chains for all cases and this largest cluster peak shifts toward right with decreasing salt due to more aggregation. The average percentage of chains, ⟨s max ⟩, belonging to the largest cluster is calculated from the broad peak in the distribution, which decreases almost linearly with increasing salt concentration as shown in the inset figure. The increase in aggregation results from the increase in average contacts per residue as shown in Figure 6(b), where the average contact for a residue is defined as the mean number of residues from the surrounding proteins within a cutoff distance 0.95 nm. The electrostatic interactions become stronger with decreasing salt concentration and the The surface mesh is constructed by alpha-shape algorithm of the OVITO visual analysis tool (probe sphere radius = 1.75σ). 43 average contact with other interchain residues shows an almost linear increase.

Figure 7(a) shows the structure factor curves for the four different salt concentrations considered. Here, each curve is suitably shifted along the vertical axis to make them distinctly visible. It can be observed that the lowest q-value peak shifts to a larger q with decreasing salt concentration and also becomes broader. Moreover, additional peaks toward the right at the higher q-values become more prominent at lower salt concentrations. These indicate more orderliness and compactness of the percolated microstructure with increasing aggregation at lower salt concentrations, which is further confirmed in Figure 7(b). We have constructed a surface mesh of the microstructure formed by the proteins using the alphashape algorithm of the OVITO visual analysis tool (probe sphere radius = 1.75σ), 43 such that the locally dense regions are represented by effective solid volumes inside the surface boundaries. The distribution of the total volume fraction of these solid volumes is shown in Figure 7(b). As seen in this figure, the effective solid volume fraction of the protein chains decreases with more amount of salt in the solution, i.e., fewer chains are part of the more compact dense regions within the percolated phase. Therefore, the system shows an overall trend of increased aggregation and compactness of the structure at lower salt concentrations.

■ SUMMARY AND CONCLUSIONS

Using a coarse-grained model of α-synuclein that accounts for its sequence, electrostatic interactions, and hydrophobicity, we have performed Langevin dynamics simulations to shed insight into the aggregation mechanism of α-synuclein. Our simulations reveal vivid details of structural evolution inside a confined volume at various protein concentrations and ionic strength. Our results complement the active pursuit reported in the literature to understand the mechanism of liquid-liquid phase separation (LLPS) in the context of formation of membraneless organelles and biomolecular condensates. Drawing an analogy with the phase behavior of associating polymers that exhibits coexistence of two phases and percolation, and mapping intrinsically disordered proteins as sticker-spacer associative polymers, the consensus in the literature is that there are two plausible mechanisms for formation of biomolecular condensates. In one scenario, LLPS occurs first, followed by percolation resulting in a gel-like phase. In the other scenarios, the protein concentration for percolation threshold is very low and the gel-like phase occurs directly without interference from LLPS. In the present study, we focus on the details of structural organization of various αsynuclein molecules at prescribed concentrations inside a fixed volume. Determination of the phase diagram for this system is relegated to future work. Our results show that at higher protein concentrations, the assembly of all chains is percolated with its structure exhibiting the earmarks of microphase separation. These microphase-separated assemblies become looser upon an increase in ionic strength. If the protein concentration is sufficiently low, then we observe only small clusters of α-synuclein molecules.

In order to further confirm that the microphase-separated percolating clusters are indeed structures in equilibrium, and not simply due to the finite size of the simulation box (∼68.4 nm, which is much smaller than a typical droplet size seen in experiments), we have performed simulations with a slab-like geometry used in literature to study LLPS in polymeric systems. 2,3 We considered a system having the dimensions 30.4 × 30.4 × 342 nm 3 with an overall protein concentration 6.44 mM and after equilibration we observe a highly concentrated region with effective concentration ∼10 mM and having microstructural organization (see Figures S9 andS10) as described above. Remarkably, the same microphase separated percolating cluster emerges in both the cubical and slab-like confinements. Even though the large scale structures at different protein concentrations presented here might suggest that the phase separation is direct gelation without interference from macroscopic LLPS, the full phase diagram needs to be computed. Nevertheless, the key conclusion from our study is that microphase separation is an integral component of assembly of α-synuclein. This important feature needs to be incorporated in constructing phase diagrams for solutions of αsynuclein as well as for other biomacromolecules.

Further analysis of our system revealed that the microstructural organization of the protein chains was mainly driven by the electrostatic interactions between the residues in Nterminal and C-terminal of the proteins. This is in tune with experimental observations of high solubility of α-synuclein at neutral pH and exposure of the N-terminus determining its aggregation propensity. 9,33 Relevance of electrostatic interactions is further disclosed when the salt concentration of the solution is varied and we find that presence of salt loosens the compactness of the microstructure as mentioned above. This can be a promising pathway to inhibit the early aggregation and ultimately amyloid formation.

Phase separation in in vivo experiments of α-synuclein occurs at much lower concentrations than in our simulations. Therefore, we believe that phase separation in solutions of αsynuclein alone might be pathological and is modulated by other PD relevant factors. 9,31 Conformational modulation of the protein due to its association with DNA has been observed 44,45 and even autogenous association with RNAs has been hypothesized. 46 How the association of these genetic biopolymers affects the phase behavior of α-synuclein is of fundamental interest and efforts in this aspect is underway. Futhermore, the liquid-like droplets ultimately mature into more solid-like amyloid fibrils, and this process is known to be influenced by the presence of lipids and other cellular components. 10 Exact conditions and interactions enabling or assisting this transition is not well understood, and insight into this is more difficult from experimental point of view. As a future step, we also intend to address this important question. Even though our simulations are centered on α-synuclein, the discovered underlying findings are a lot more general in the context of biomolecular condensates.