Postsynaptic protein assembly in three- and two-dimensions studied by mesoscopic simulations

Recently, cellular biomolecular condensates formed via phase separation have received considerable attention. While they can be formed either in cytosol (denoted as 3D) or beneath the membrane (2D), the underlying difference between the two has not been well clarified. To compare the phase behaviors in 3D and 2D, postsynaptic density (PSD) serves as a model system. PSD is a protein condensate located under the postsynaptic membrane that influences the localization of glutamate receptors and thus contributes to synaptic plasticity. Recent in vitro studies have revealed the formation of droplets of various soluble PSD proteins via liquid-liquid phase separation. However, it is unclear how these protein condensates are formed beneath the membrane and how they specifically affect the localization of glutamate receptors in the membrane. In this study, focusing on the mixture of a glutamate receptor complex, AMPAR-TARP, and a ubiquitous scaffolding protein, PSD-95, we constructed a mesoscopic model of protein-domain interactions in PSD and performed comparative molecular simulations. The results showed a sharp contrast in the phase behaviors of protein assemblies in 3D and those under the membrane (2D). A mixture of a soluble variant of the AMPAR-TARP complex and PSD-95 in the 3D system resulted in a phase-separated condensate, which was consistent with the experimental results. However, with identical domain interactions, AMPAR-TARP embedded in the membrane formed clusters with PSD-95, but did not form a stable separated phase. Thus, the cluster formation behaviors of PSD proteins in the 3D and 2D systems were distinct. The current study suggests that, more generally, stable phase separation can be more difficult to achieve in and beneath the membrane than in 3D systems. SIGNIFICANCE Synaptic plasticity is a key factor in memory and learning. Upon learning, protein condensates that form beneath the postsynaptic membrane are known to change their nature. Recent studies have suggested that condensate formation is related to liquid-liquid phase separation based on in vitro experiments of soluble parts. However, the phase behavior can be strongly dependent on physical dimensions. The mechanism by which condensate grows beneath the membrane is not well characterized. Taking advantage of the ease of systematic comparison using computer simulations, we investigated the phase behaviors of postsynaptic protein assemblies in 3D and 2D systems. The results revealed that even when a 3D system exhibited clear phase separation, the corresponding 2D system did not exhibit it stably.


INTRODUCTION
Recently, biomolecular condensates in cells have received considerable attention (1).
Biomolecular condensates are heterogeneous assemblies that are composed of many proteins and RNAs. As these condensates are considered to have cellular functions but are not surrounded by membranes, they are often called membrane-less organelles. Classic examples include Cajal bodies, nuclear speckles, stress granules, and postsynaptic density (PSD), whereas recent studies have revealed many more examples, such as signalosomes (2). These condensates are formed mainly by two types of interactions: interactions mediated by disordered regions with less sequence specificity, and specific and stoichiometric interactions between protein domains (3)(4)(5)(6)(7). Notably, these condensates can form either in solutions, as in the case of Cajal bodies (the 3D system), or under the membrane, as in the case of PSD (5,8) (the 2D system). Biomolecular condensates often form via liquid-liquid phase separation (LLPS). It is well-known in statistical physics that phase transition behaviors are strongly dependent on the underlying physical dimension (9); 1D systems cannot show any phase separation, 2D systems are marginal, and 3D systems tend to exhibit sharp phase transitions. This trend may sound odd because LLPS in biology has often been reported in 2D systems. While 3D and 2D LLPS behaviors can qualitatively be different, the possible differences have not been well explored, to the best of our knowledge.
To compare the phase behaviors in 3D and 2D, PSD proteins can serve as a model system. PSD is a protein condensate beneath the postsynaptic membrane that was first observed by electron microscopy at the distal tip of the spine head (10). In hippocampal CA1 pyramidal cells, PSD has an proposed domain-resolution mesoscopic model makes comparative simulations of the LLPS dynamics tractable (30,31). In this model, specific interactions between protein domains are treated as virtual chemical bonds formed by reactions, which ensure one-to-one stoichiometric binding within the mesoscopic representation. The key parameters in the model were calibrated using a 3D system based on in vitro experimental data. We analyzed the protein domain interaction networks formed in the 3D and 2D systems, and found significantly reduced network edges in the 2D system. While the 3D system exhibited clear phase separation, scaling analysis of the 2D system indicated that the mixture of AMPARs, TARPs, and PSD-95 beneath the membrane formed finite-sized clusters, but not the separated phase. In general, the results suggest that the phase separation behaviors are markedly different between the 3D and 2D systems. The droplet formation in a soluble 3D system does not directly prove phase separation in the corresponding 2D system.

Mesoscopic model: Molecular representation
Based on previous experimental assays (28,29), we simulated the mixture of a glutamate receptor AMPAR (or its soluble counterpart), a regulatory protein TARP, and a scaffolding protein, PSD-95, using a mesoscopic model (30,32,33). In the mesoscopic model, each globular domain was represented by a single spherical particle. The neighboring domains are connected by interactions derived from the Gaussian polymer model, which represents flexible linkers.
In this study, we considered two simulation systems: 1) a system containing soluble variants of AMPAR and TARP, together with PSD-95s. AMPAR and TARP were replaced with DsRed2 (PDB ID: 1ZGO, N = 936 residues) and the soluble C-terminal part of TARP (TARPc), respectively, based on the experiment (29) (denoted as the 3D system) (Fig.1A top); and 2) the system that contained AMPAR and TARP restrained to the membrane surface, in addition to palmitoylated PSD-95 (denoted as the 2D system) (Fig.1B top). Thus, in addition to the underlying geometry, the two systems have different molecular sizes and diffusion coefficients. However, given the mesoscopic nature of the model, these differences would not severely alter the results.
Although AMPAR forms a tetrameric structure constituted by GluA1-4 subunits, in this study, we regarded it as a single sphere, assuming that its tetrameric structure was maintained throughout the simulation. The AMPAR is assumed to be composed of a homo-tetramer of the GluA2 subunit (PDB ID: 3kg2, N = 3292 residues). TARP is a regulatory protein that binds to the GluA subunit on a one-to-one basis, thereby regulating AMPAR protein levels, synaptic localization, and channel activity (34). As a representative of TARP proteins, we chose the TARP γ 2 subunit (N = 211 residues), also known as "stargazin" since it has been studied as the most well-known TARP subtype in numerous preceding researches. We modeled the TARP (or TARPc) protein as a single spherical particle. Throughout the simulation, we assumed that four TARP particles were always connected to one AMPAR (or DsRed2) and never dissociated from AMPAR. It can thus be denoted as DsRed2(TARPc) 4 and AMPAR(TARP) 4 for the 3D and 2D systems, respectively. PSD-95 contains five globular domains: three PDZ domains (PDZ1, PDZ2, and PDZ3, all N = 86 residues), one SH3 domain (N = 70 residues), and one GK domain (N = 165 residues) (28,35). In addition, an additional small particle representing the first five amino acids of PSD-95 with two cysteine residues, which are palmitoylated in the case of 2D system, is added to the N-terminus (N = 5 residues) and tethered to the postsynaptic membrane (36). We represent a PSD-95 molecule as six spherical particles connected by a Gaussian chain. In the 3D system, a single sphere representing an AMPA receptor in the 2D system was substituted by a soluble fluorescent protein called DsRed2, which also has a tetrameric structure (29). The N-terminal particle of PSD-95 was also replaced with a sphere, with no restraint on the membrane.
The radii ܴ (in nm) and diffusion coefficients ‫ܦ‬ of globular domains were set using the formula below as a function of the number of amino acids ܰ in the domain (37): is the hydration radius, ݇ is the Boltzmann constant, ܶ is the temperature, and η is viscosity. Since ܴ is determined by the number of amino acid residues, molecular size of AMPAR in the 2D system is larger than that of DsRed2 in the 3D system for the representation of a tetrameric protein with 4 TARPs. Viscosity in the 3D system and that of the domains in the cytosol in the 2D system were determined by water viscosity at 300 K (0.89 cP). The viscosity of the membrane-embedded domains in the 2D system was set to be 100 times that of the cytosol, since the viscosity of the bilayer core in the DPPC vesicle was estimated to be 80-100 cP (38).

Mesoscopic model: Potentials
Based on the molecular representation in each system, the molecular system is subject to the total potential energy function as follows: and interactions derived from the Gaussian polymer chain . We set the bond between AMPAR and TARP as a harmonic bond. The spring coefficient of the harmonic bonds was set to 10 kJ/mol/nm 2 . On the other hand, protein linkers in PSD-95 are represented by a Gaussian polymer chain, which is easy to extend and shrink, reflecting the flexibility of the protein linker. Compared to that of the harmonic bond, which is set relatively high to maintain a stable tetrameric structure, the spring constant of the Gaussian chain was determined by the number of amino acids (ܰ) contained in the linker and takes the value from k ~ 0.3 to 3 in the current case.
was set to 10 kJ/mol/nm 2 , which is the same as that for other harmonic bonds. Here, we call the interaction defined as a virtual bond as a "specific interaction." By appropriately defining virtual reaction systems, we can rigorously exclude the binding of PDZ domains to TARPs that are already bound to other PDZ domains, guaranteeing a one-to-one specific interaction. Although the C-terminus of TARP, which contains a PDZ binding motif, is said to interact mainly with PDZ1 or PDZ2 on the postsynaptic membrane (28), we assumed that all PDZ domains have the same binding ability to TARPc in this model, as the binding ability to PDZ3 has also been confirmed by artificially extending the C-terminal tail of TARP (35).
In addition to the virtual bond created by the virtual reaction, we implemented the nonspecific interaction ܸ ௦ which acts between any two protein domains, in contrast to the specific interaction. This interaction consists of a weak inter-domain force and an excluded volume term between the two particles. The potential function is obtained by connecting the three harmonic functions as follows (Fig.1D) , was tuned based on the experimental data, as described below. Note is set to zero, a repulsive potential is applied when the intermolecular distance is smaller than was set to 2.5 nm. To achieve this, we introduced an additional potential that implicitly represents the membrane. We used transmembrane domains (Fig.1E, right). The time step of the simulation was 1.0 ns for the 2D system and 0.25 ns for the 3D system and the parameter tuning simulation. With time steps larger than these, we occasionally observed breakage of bonds between AMPAR and TARP. Given the simple Langevin dynamics used, the time scale in the software does not necessarily correspond precisely to actual physical time.
Simulation Protocol: Slab simulation in the 3D system interactions. The same simulation protocol was used to simulate LLPS with both specific and non-specific interactions.
As described above, the specific interaction between TARPc and PDZ is represented by a virtual reaction such that close molecules are tied together by virtual bonds. Here, we call the groups of molecules linked by virtual bonds as "clusters." Simulation Protocol: Parameter tuning in the 3D system We tuned the parameters for the specific and nonspecific interactions (ܸ ) based on experimental data. In the specific interaction, we fixed ݇ and the threshold distance and altered ݇ to control the specific interaction strength, . In nonspecific interactions, we tuned the depth of the potential well, are the cumulative durations of the unbound state, bound state, and total simulation time, respectively.
ܸ is the volume of the simulation box.
ܰ is Avogadro's constant. In this tuning step, we assumed that two molecules are in "bound state" when the distance between TARPc and any PDZ domain is smaller than 1 .

ԛ ‫ݎ‬
, which is the reference length at which two particles begin to interact. Conversely, "unbound state" is when the distances between TARPc and every PDZ domain are larger than the interaction length. This definition of bound/unbound state does not depend on the formation of a virtual bond created by specific interactions.
The second assay measured the concentration of the dilute phase when the system exhibited phase separation during the slab simulation. Experimentally, the critical concentration for phase separation, which corresponds to the concentration of the dilute phase in the phase-separated system, was ‫ػ‬ 1 0 µ M . We note that DsRed2(TARPc) 4 and PSD-95 have four and three interacting domains that interact one-to-one, respectively, and that the simulation box contains 200 molecules of each.
This indicates that there is an excess of DsRed2(TARPc) 4 as a client for specific interactions in this system. Therefore, we used the concentration of PSD-95 in the dilute phase as a reference for the dilute phase concentration. Adopting the slab simulation protocol described above, we repeated the molecular simulations with different parameter sets that provided appropriate In the MSD calculation, we first separated the slab box into condensed and dilute phases as described above. Then, in each phase, we obtained the MSD of the molecules as a function of time.
Even though we used the periodic boundary condition, we could obtain the absolute distance of diffusion by tracking the motion. As in the FRAP experiment, we took an average of 1000 snapshots over the reference time point to obtain smooth MSD lines. The diffusion coefficients of the molecule in each phase were estimated using the value calculated from the slope of the MSD.
Simulation Protocol: Simulation in the 2D system

RESULTS AND DISCUSSIONS
Model validation and phase separation in the 3D system We first examined our mesoscopic modeling for the mixture of PSD-95 and DsRed2(TARPc) 4 , a soluble counterpart of AMPARs bound to the four C-terminal cytoplasmic parts of TARPs, based on available experimental data. This system has recently been shown to exhibit LLPS in in vitro assays.
The C-terminus of stargazin (TARP We adopted the slab simulation box, which is known to be advantageous for realizing phase separation in MD simulations (42,43). With more standard cubic boxes, protein condensates would take near-spherical droplets which are entirely surrounded by interface. The free energy cost from the interface results in large finite-size artifacts. The slab shape with the periodic boundary condition enables the separation of the two phases by the interfaces of small areas, which alleviates the finite-size effects (42,43). We began the simulations using homogeneous and random configurations.
We first tested modeling solely by the specific one-to-one interaction between each of the PDZ domains of PSD-95 and TARPc. We set the nonspecific interaction to zero, i.e., formed small clusters, which gradually coalesced into larger clusters (a final snapshot in the second panel of Fig.2A). At the end of the 20 ms (8 × 10 7 MD time step) simulations, we found no molecules in the dilute phase. We note that the concentration of the dilute phase corresponds to the critical concentration for phase separation (44), which indicates that the critical concentration of LLPS in this setup is much lower than 1 μ M. This is inconsistent with the experimental results.
We then explored the two-dimensional parameter space kJ/mol, the system did not maintain the phase separation resulting in uniform distribution (Fig.2B) Fig.3A and movie in the supplemental data, Movie S1). Figure 3B shows the time course of the number of virtual bonds, i.e., the one-to-one binding between TARPc and PDZ (blue) and the number of DsRed2 (red) in the largest cluster. Since the time course of the cluster size has reached a plateau, we considered that the simulation reached the equilibrium at ~40 ms. The increase in virtual bonds corresponds to the progress of domain-wise protein network formation, whereas the largest cluster size represents the progress of condensed phase formation. The state of the virtual bonds converged to equilibrium much faster than the phase behavior. As the fusion of two droplets doubles the number of particles in the largest cluster, we observed that the largest cluster grew step-by-step during the simulation. The coefficient of variation of the largest cluster size (Fig.S1B) shows that once a droplet is formed, it maintains the composition stably with little fluctuation.
Notably, after the coalescence of two spherical droplets at early stages, the double-sized droplet rapidly assumed a nearly spherical shape. This is a hallmark of droplet fluidity.
To further assess whether the formed condensate has fluidity, we conducted fluorescence Cluster formation on and beneath the membrane (the 2D system).
To address the difference between the assembly of the soluble counterpart of AMPAR, i.e., DsRed2, in in vitro experiment and that of AMPAR on the postsynaptic membrane, we performed simulations in a 2D system. In contrast to DsRed2(TARPc) 4 , AMPAR and TARP particles were restrained to the membrane. In addition, the N-terminal particle of PSD-95 was assumed to represent the palmitoylated group and was thus restrained to the membrane. These membrane-embedded particles have smaller diffusion coefficients in the membrane than their soluble counterparts in the solution.
For interactions between TARP and the PDZ domain of PSD-95, we used the same parameter set as that used in the 3D system.
We first performed molecular simulations with 50 PSD-95 and 50 AMPAR(TARP) 4  adult rat cerebellum (51,52), respectively. A dense condition was used to investigate the limiting case. Figure 5A shows the final snapshots of the simulation from a random and homogeneous configuration in the three densities from two directions, viewpoints top and side in Fig.1B (bottom).
The results showed that PSD-95s assembled AMPARs to form clusters on the membrane at any areal density. Viewpoint side shows that the formation of AMPAR clusters (red particles) was co-localized with the PSD-95 molecules beneath the clusters. The temporal changes in the largest AMPARs clusters under three different areal density conditions (Fig.5B) showed that the lower the density, the AMPARs belonging to large clusters, single or small groups of AMPARs were also observed on the membrane. The features of cluster formation in the 2D system, such as stepwise growth of the cluster and slower domain-wise binding, were shared with all three areal density cases.

Receptor-scaffold network does not exhibit the phase separation on the membrane
Although both the mixture of DsRed2(TARPc) 4 and PSD-95 in the 3D system and that of AMPAR(TARP) 4 and PSD-95 in the 2D system formed clusters, we also observed that the stability of the large clusters differed between the two systems. Large clusters formed on and beneath the membrane were fragile, and breakage was often observed. On the other hand, we observed virtually no breakage of large clusters formed in the 3D system. In addition, the observed cluster sizes in the 2D system differed among the areal densities: a large cluster of receptors constituted by almost all the receptors in the system was observed in the dense condition but not in other conditions. While phase separation was clearly visible in the 3D system, it is unclear whether the cluster formation of receptors and scaffold proteins in the 2D system can be regarded as phase separation.
To address the phase behaviors more rigorously, we examined the scaling behavior both in the 3D and 2D systems. We note that for a cluster (droplet) formed in a multicomponent system to be regarded as a "phase," the size of the phase should grow linearly with the system size in a set of self-similar systems (53,54). Thus, we examined the self-similarity of droplets in a series of 3D and 2D systems of different sizes with identical densities (In the 2D system, we employed the dense condition). Both in 3D and 2D systems, we repeated simulations for the mixture of AMPAR-(TARP) 4 and PSD-95 for four different numbers of AMPARs and PSD-95s at 25, 50, 100, and 200, and measured the number of AMPARs in the largest cluster. We note that, in the 3D system, while AMPAR is all substituted by DsRed2, we denoted it as AMPAR for simplicity.
To check the convergence of simulations, which was not easy particularly for larger systems, we employed a pair of simulations for each system size: a run from a random and homogeneous configuration (termed "forward") and another run from a perfectly phase-separated configuration cluster size after merging the two trajectories (the right panels in Fig.6A) and the evolutions of coefficient of variation (Fig.S2E) show broad distributions in the 2D system, indicating large intrinsic fluctuations. In the 2D system, snapshots in the final simulation step (Fig.6B right) suggest that the molecules form essentially one cluster for the case of 25 AMPARs but have several clusters for the case of 200 AMPARs. This is in sharp contrast to the geometrically similar growth of one large condensate in the 3D system.
Finally, we plotted the average size of the largest AMPAR cluster as a function of the number of AMPARs in the system (Fig.6C). The results indicate that although the largest AMPAR cluster size increases in both systems, the increase is not proportional to the number of AMPARs in the 2D system, contrary to the linear growth in the 3D system. The largest AMPAR cluster size in the 2D system seemed to saturate below ~100. Sublinear growth and plausible saturation of the cluster size are hallmarks of the lack of phase separation. To assess the finite-size effect of the simulation box, we also performed a 2D simulation in a 2D-variant of the slab box finding that the scaling of the largest AMPAR cluster size is not significantly affected by the different box setup (Fig.S3). We also confirmed that there is little effect on the scaling of the largest AMPAR cluster size even under the planar constraints to the geometry of four TARPs, by introducing additional bonds and dihedral potentials in AMPAR(TARP) 4 (Fig.S4).
We note that, in the 2D system, the areal density under this condition is probably larger than that of AMPARs on the actual postsynaptic membrane, as mentioned above. Even under such over-dense conditions, stable phase separation was not observed. With physiological AMPAR density, phase separation of the mixture of AMPAR(TARP) 4 s and PSD-95s is less likely to occur.
Therefore, we suggest that any AMPAR density found in the postsynaptic membrane, such as in the hippocampus or cerebellum, is unlikely to cause stable phase separations on and beneath the membrane with only the binary mixture of AMPAR(TARP) 4 s and PSD-95s. Indeed, in an experiment in which receptor-like proteins were tethered on one side of the lipid bilayer to mimic PSD diffusion beneath the membrane, removal of even one of the four scaffolding proteins (PSD-95, GKAP, Shank3, Homer3) was found to completely prevent PSD cluster formation, supporting our simulation results (8).

Differences in network structure between the 3D and 2D systems
Finally, we investigated the multivalent interactions between AMPAR (or its soluble counterpart, DsRed2. Henceforth, we denote it AMPAR for simplicity), and PSD-95 via TARPs and PDZ domains leads to differences in the network structure between the 3D and 2D systems. As the specific interactions set up throughout this study are treated as virtual association and dissociation reactions, it is possible to specify which TARP specifically interacts with the PDZ domain. Using this property, we examine the distribution of virtual bonds.
We counted the average number of PDZ domains bound per AMPAR from the final 1000 snapshots of the simulation, with 200 of both molecules added to the simulation box. There was a somewhat interesting difference in distribution between the two systems: the results showed that the four TARPs per AMPAR had the highest population in the 3D system, whereas three PDZ domains per AMPAR were found to be the most common in the 2D system (Fig.S5A). This suggests that the protein-domain network in the 3D system exhibits a dense and complicated entangled structure beyond the one-to-one binding of AMPAR complexes and PSD-95s.
To maximize the effect of multivalent interactions, which are the main driving force of phase separation, TARPs were assumed to bind to all four AMPAR subunits in the current work. However, since AMPARs have functional regulatory subunits other than TARPs that competitively bind to the same structural positions as TARPs, this assumption of AMPAR with four TARPs may be unlikely to occur at real synapses (55)(56)(57). The results of the virtual bond network analysis imply that even under the condition where the maximum number of TARP (four) is available as an interactive particle, three TARPs are sufficient for AMPARs to form collective structures on the actual postsynaptic membrane, owing to the spatial constraint imposed by the membrane.
The distribution of the number of TARPs bound per PSD-95 also showed slight differences between the two systems: the number of TARPs bound per PSD-95 was almost three in the 3D system, while two were also observed to no small extent in the 2D system (Fig.S5B). Generally, PSD-95 is more readily incorporated into the condensed phase than AMPAR in both systems as there are four TARPs per AMPAR and three PDZ domains per molecule of PSD-95. The setting of the number of molecules may indicate that PSD-95 is relatively more likely to bind more TARPs than the tendency of PDZs to be bound by AMPAR (Fig.S5A), which is consistent with the results for the number of bound TARPs in 3D. However, a surplus of PDZ domains that can interact with TARP was identified in the 2D system despite TARP outnumbering PDZ in the system. Similar to the AMPAR results, this result may be explained by the fact that the number of TARPs that each PSD-95 can access is spatially confined by the geometry of the postsynaptic membrane.
The multiplicity of virtual bonds between a pair of AMPAR and PSD-95 also revealed dimensional differences. The specific interaction between AMPAR and PSD-95 can form up to three virtual bonds (Fig.S5C). The valences of the bonds observed from all pairs of bound AMPAR and PSD-95 pairs were nearly monovalent in both systems, indicating that multi-bonding within the two molecules was relatively rare. However, it was found that the bond valency tended to be higher in the 2D system than in the 3D system (Fig.S5C).
The characteristics of the virtual bond networks in the 3D/2D systems shown in Figs.S4A, B, and C become more prominent when compared with the graphical diagrams ( Fig.S5D and S5E). Figure S4D and S4E illustrate the local bond networks from the final snapshots in the 3D and 2D systems, respectively. These diagrams were drawn starting from AMPAR in the center and tracing 11 bonds (regardless of bond type). For particles shown with a single star, the AMPAR is closest to the center of gravity of the condensed phase in 3D or the largest cluster in 2D. The graphs show that almost all PDZ domains (yellow) are occupied by specific bonds, whereas some TARPs are not used for binding in either system. Notably, the network structure is more developed and densely formed in a 3D system than in a 2D system. The network diagram does not clearly show the difference in the multiplicity of virtual bonds as the formation of multivalent bonds seems to be a relatively rare in both the systems (Fig.S5C). However, the difference in the development of a network is presumed to be correlated with binding multiplicity: in a 3D system, a single molecule, such as AMPAR or PSD-95, can interact with many other molecules. In contrast, the membrane in 2D systems restricts the existing area such that the molecules in and below the membrane can only interact with each other, thus preventing the formation of a dense network. This difference in network structure is hypothesized to be due to the special constraints imposed by the membrane, and it may be a reason why phase separation is less likely to occur in the 2D system than in the 3D system. PSD is a protein condensate formed beneath the postsynaptic membrane and is correlated with signal transmission and synaptic plasticity. This molecular architecture beneath the postsynaptic membrane has been observed using electron microscopy, suggesting that the condensate is stable (10,16). However, the current simulation results suggest that AMPAR(TARP) 4 and PSD-95 alone form transient small clusters, but are not sufficient to form a stable phase immediately under the membrane. These are not contradictory as PSD is known to be organized in a hierarchical manner with other scaffolding proteins, such as SAPAP, Shank, and Homer, which are also prone to form protein condensates in in vitro 3D systems (8,58,59), other than PSD-95, which is generally located in the upper hierarchy of PSD. At the bottom of the PSD hierarchy, the actin cytoskeleton further contributes to retaining the protein network and anchors the condensate to the cell body. With this hierarchy, protein networks may physically be between 2D and 3D systems, which could enhance the stability of the condensate. Incorporation of these additional scaffolding and cytoskeleton proteins into mesoscopic computational modeling would be future challenge.

CONCLUSIONS
More generally, we suggest that qualitative aspect of the current results is applicable to other protein assemblies with multidomain interactions, since the mesoscale model we employed only depends on a limited number of parameters and is less protein-specific. Further research may explore a more general consideration of the extent to which domain multi-valency and protein-protein affinity affect phase separation.
In general, the phase transition depends on physical dimensions (9,60,61

DECLARATION OF INTERESTS
The authors declare no competing interests. Figure S1. Convergence examination. Figure S2. Preparation of the initial configuration in the reverse simulation in Fig.6 and the coefficients of variation in two systems. Figure S3. Simulation in a 2D system with a slab box. Figure S4. Size scaling analysis with planar geometry.      C) The mean and the standard error of the largest AMPAR cluster size for the three areal densities.