Quantitative analysis of membrane-mediated interactions and aggregation of �exible peripheral proteins orchestrating large-scale membrane remodeling

Shaping and remodeling of biomembranes is essential for cellular trafﬁcking, with membrane-binding peripheral proteins playing the key role in it. Signiﬁcant membrane remodeling as in endo-and exocytosis is often due to clusters or aggregates of many proteins whose interactions may be direct or mediated via the membrane. While computer simulation could be an important tool to disentangle these interactions and understand what drives cooperative protein interactions in membrane remodeling, this has so far been extremely challenging: protein-membrane systems involve time-and lengthscales that make detailed atomistic simulations impractical, while most coarse-grained models lack the degree of detail needed to resolve the dynamics and physical effect of protein and membrane ﬂexibility. Here, we develop a coarse-grained model of the bilayer membrane bestrewed with rotationally-symmetric ﬂexible membrane-bound proteins. We show how this model can be parameterized based on local curvatures, protein ﬂexibility, and the in-plane dynamics of proteins. We measure the effective interaction potential for the membrane-mediated interactions between peripheral proteins. Furthermore, we investigate the kinetics, equilibrium distributions, and the free energy landscape governing the formation and break-up of protein clusters on the surface of the membrane. We demonstrate how the ﬂexibility of the protein plays a deciding role in highly selective macroscopic aggregation behavior. Finally, we present large-scale simulations of membrane tubulation, and discuss the sequence of events and the stability of intermediates.

Membrane-bound peripheral proteins are key players in a wide range of biological functions, especially when they are capable of sensing/inducing membrane curvature.These proteins are often found in regions of the plasma membrane that undergo significant conformational change in processes such as tubulation, exo/endocytosis, and fission/fusion.The same also holds for pathogens such as Shiga and Cholera toxins [1][2][3].It was established early on that the internalization of Cholera toxin from cell membrane into endoplasmic reticulum can progress with or without relying on clathrin and caveolae, i.e. active endocytosis [4].The pentagonal subunit B of Shiga toxin (STxB) was also shown to induce invaginations, and long tubules, in the absence of active cellular machinery [2].This rather elegant entry mechanism has even inspired using STxB as a vehicle for drug delivery and anticancer therapy purposes [5,6].Considering the rather large bending stiffness of bilayer membranes (∼20 kT), macroscopic membrane remodeling involves significant energy expenditure that is beyond the action of any single membrane-binding protein, thus necessitating the cooperative action of an assembly of peripheral proteins [7].Independent of binding and curvature-generation mechanisms [8][9][10][11], the concerted action of a collection of these proteins can lead to macroscopic membrane deformations [12].
Due to limitations of simultaneously resolving time and space at the nanoscale with experimental appraoches, physics-based computer models are a key technology to gain insights into structural and spatiotemporal mechanisms.However, physics-based models of the peripheral protein clusters are challenged by multiscale nature of the problem: on one hand, a simulation model must provide enough detail to resolve the relevant molecular effects, on the other hand, it needs to cover time-and lengthscales not simultaneously achievable by atomistic models.
For example, lipid reorganization has been indicated when these proteins are in close proximity [13].Yet, macroscopic dynamics of protein clusters implies the existence of long-range membrane-mediated interactions (see [14] and references within).Studying the nature of these interactions has a long history in membrane biophysics [15][16][17][18][19][20][21].Simulation tools have proved essential in the study of these interactions and the resulting organization [22][23][24][25][26]. Non-symmetric proteins such as crescent-shaped BAR-domains experience strong membrane-mediated interactions resulting in salient aggregation and structuring on the membrane [12,[27][28][29].However, the interactions and cooperativity of rotationally symmetric proteins, such as STxB, rely on the much more subtle fluctuation-induced forces [14][15][16][17][18].These so-called thermal Casimir interactions bring proteins that suppress membrane fluctuations together to minimize their adverse entropic effect.Despite several studies, and a multitude of hypotheses, reliable quantitative results that offer insight into the cooperative action of these proteins, in well-defined pathways resulting in membrane remodeling, are scarce.Also, whilst the kinetics of membrane and peripheral proteins are of utmost importance in remodeling processes [30], developing a model that incorporates in-plane kinetics of the proteins on the same grounds as the large-scale kinetics of the membrane and solvent is still challenging.
Here we develop a coarse-grained model of flexible membrane-bound proteins interacting with dynamic membranes.The model is fully particle-based so as to facilitate the integration with particle-based simulations of cellular kinetics [31,32], and particularly with interacting-particle reaction dynamics (iPRD) [33][34][35][36].Our model builds on a previous coarse-grained particlebased membrane model mimicking the mechanics of the bilayer, while coupling to the hydrodynamics of the solvent to reproduce realistic large-scale equilibrium and non-equilibrium kinetics [37][38][39].This model is extended to incorporate cooperative membrane-mediated interactions as well as the dynamics of protein aggregation, in order to arrive at a realistic model of the dynamic membrane-protein interplay.Our model can be parameterized solely with experimental observables, such as the bending rigidity and the solvent viscosity.Using large-scale simulations, we extract potential functions for indirect membrane-mediated interactions between peripheral proteins.We also investigate the thermodynamics and kinetics of protein aggregation on the membrane and demonstrate opposing effects of protein stiffness.We suggest, somewhat contrary to previous studies [26], that optimal aggregation does not occur for stiffest proteins/strongest membrane binding but at an intermediate flexibility.We also investigate the free energy landscape of tubulation process as a result of such aggregates, to reveal very large barriers and highly stable intermediates, offering elaboration of the hypotheses regarding the membrane remodeling following spontaneous protein aggregation [40].

Particle-based membrane model
For the simulations presented here, we extend the particle-based membrane model of ref. [37], composed of particle-dimers with nearest-neighbor bonded interactions.We have chosen a lattice parameter of 6.5 nm, which is comparable with the radius of the STxB protein, to represent a membrane that allows tightly packed clusters of proteins.Membrane particles can be tagged to represent where peripheral proteins are bound (Fig. 1).All effects corresponding to the presence of a bound protein are incorporated in a locally masked force field.We utilize the combined energy density of membrane and the protein (details in Methods section "Elastic model of membrane and peripheral proteins") in conjunction with parameter-space optimization (details in Supplementary Information) to obtain force field parameters for the unmasked (membrane) and masked (membrane + protein) regions (inset plots in Fig. 1).

Effective membrane curvatures
We can obtain an initial estimate of the local membrane curvature using the distribution of energy density in membrane + protein regions, f eff (Fig. 1 inset plots).This value depends on the preferred curvature of the protein as well as the stiffness attributed to it (Methods section "Elastic model of membrane and peripheral proteins").The actual curvature observed in simulations would deviate because: (i) the parameterization scheme is only applied considering the nearest neighborhood of each particle, and (ii) in order to observe detailed balance, the local masking of the force field is achieved through Monte Carlo moves (Methods section "Simulations").
We have simulated ∼325×325 µm 2 membrane patches containing 10 % surface concentration of bound proteins with stiffness Y p = 50, 100 and 200 MPa (details in Methods section) to measure actual curvatures.As expected, stiffer peripheral proteins lead to more pronounced membrane curvatures (Fig. 2) with the energy-based estimate of curvature providing an upper bound (Fig. 2, arrows).The wide range of curvatures facilitates modeling a variety of proteins (Fig. 2).We can compare these curvatures to the reported results from molecular dynamics simulation of STxB bound to Gb3 lipids, which yields inward membrane curvatures of 0.034 ± 0.004 nm −1 for saturated and 0.035 ± 0.003 nm −1 for unsaturated acyl chains on the lipids [41].Also, using curvature sorting of Cholera toxin subunit B (CTxB), the induced membrane curvature has been experimentally determined to be 0.055 ± 0.012 nm −1 [42].All these values are within the range offered by proteins studied here (Fig. 2).

Kinetics of membrane and proteins
In order to predict realistic timescales and mechanisms, the model must faithfully capture the kinetics for the membrane and membrane-bound proteins.We have previously done that for the kinetics of out-of-plane membrane fluctuations with the current model [39].Focusing on the in-plane kinetics, here we investigate the in-plane diffusion of protein-bound particles (Fig. 3a) and subsequently calibrate it to the experimental range of values.In our model, two control parameters, respectively the frequency of bond-flipping moves and the inplane mobility of particles, control the diffusion (Methods section and [37,39]).
We observe the mean-squared displacement (MSD) to follow the familiar three-step regime [43]: fast diffusion for short time intervals, followed by subdiffusive behavior, signifying the onset of crowding effects, and finally, the asymptotic normal diffusion for large lag times (Fig. 3a).
Our model reproduces the expected 4D p t behavior for the asymptotic diffusion regime (Fig. 3a).We have been Figure 1: Schematic of the model: mesoscopic membrane model with some particles "tagged" as supporting peripheral membrane-bound proteins.Where a protein is bound to the membrane, a locally modified force field is applied (masking of the force field).Force field parameters are optimized to reproduce the Helfrich energy density for regular particles, plus an additional strain energy corresponding to the deformation of peripheral proteins for tagged particles (details in Methods section "Elastic model of membrane and peripheral proteins" and Supplementary Information).The two plots compare the energy densities of the theoretical models and the effective energy densities resulting from particle-particle interactions with the optimized parameters for the given range of membrane curvatures.In the lower plot, results are given where the Young modulus of Y p = 100 MPa, and a preferred curvature of c 0 = 0.08 nm −1 are attributed to the peripheral proteins.
able to obtain protein diffusion coefficients comparable with the experimental measurements (Fig. 3a).Using fluorescence recovery after photobleaching (FRAP) assays, Tian and Baumgart measured a diffusion coefficient of 0.35 ± 0.09 µm 2 s −1 for CTxB proteins [42].Also, Day and Kenworthy similarly measured the diffusion coefficient of STxB proteins as ∼0.5 µm 2 s −1 [44].It is noteworthy that multivalent glycolipid-binding proteins such as CTxB or STxB diffuse extremely slowly compared to their lipid-anchored counterparts [45].The in-plane diffusion coefficients show a slight dependence on the protein stiffness, with Y p =100 MPa proteins diffusing more slowly, while the mentioned control parameters are similar for all cases (Fig. 3a).
We obtain the macroscopic surface viscosity of the membrane, based on fluctuations of the in-plane shear stress and the corresponding Green-Kubo relation [38,39,46,47].Irrespective of the protein stiffness, timecorrelation functions of in-plane shear stress decay very rapidly and do not show any long-term memory effects (Fig. 3b).Surface viscosity, especially considering the error margins, only weakly depends on the protein stiffness, with slower kinetics observed again for Y p =100 MPa proteins (Fig. 3b).The experimental value for the surface viscosity of the membrane is in the 10 −7 to 10 −6 sp range [48][49][50][51].While molecular dynamics simulations with coarse-grained force fields usually result in much lower viscosities in the 10 −10 to 10 −8 sp range [52,53].Despite being highly coarse-grained, our model produces values in much better agreement with experiments.The timescales obtained for diffusionlimited processes will thus be reliable.

Membrane-mediated interactions
Contrary to previous models with similar scale and resolution (see for example [41]), we have not included any direct interactions between peripheral proteins.We are thus able to observe and measure the emergent membrane-mediated interactions.To this end, we employ the statistical mechanics of liquid mixtures to estimate the effective interactions between pairs of the constituents.We consider a two-dimensional fluid, consisting of either membrane (m) or protein-bound (p) particles, and respectively denote the effective pairwise interactions as u mm , u mp , and u pp .These effective interactions are estimated based on pair-correlation functions (details in Methods section "Statistical mechanics of liquid mixtures").The interaction potential u pp also contains contributions from the underlying membrane particles.To correct for this, we have simply calculated the net membrane mediated interaction as u pp −u mm + 1 2 u mp (Fig. 4a).
We consider two prominent continuum-based models for comparison: (i) repulsive interactions between symmetrically curved regions, either from peripheral proteins or from stiff inclusions in the membrane [7].These interactions stem from the membrane elastic energy (Eq.( 4)) acquiring its minimum only when the two curved regions are infinitely apart.For a pair of particles with the separation r, this interaction is approximated in the leading power by [15,18], where θ p and r p are the contact angle and the radius of the curved region and κ is the bending rigidity of the membrane (Eq.( 4)).(ii) attractive entropic forces emerging when the binding is stiff enough to locally suppress membrane fluctuations (the so-called thermal  Casimir effect).An approximate expression for the potential of this attractive interaction is [15,18,20],   1) and ( 2)).Contact angles (θ p ) are calculated based on the measured mean curvatures given in Fig. 2. (b) Effect of the protein stiffness on the shape and strength of the interaction potential.Well-depth, ∆u w , and barrier height, ∆u b , are measured according to the inset schematic shown in (a).The red and gray shaded regions in (a) are where the local minima and maxima of potential energies have been sampled.
Membrane-mediated interaction potentials that we obtained from the simulations show significant attraction at close range (Figs. 4a).We also observe a long-range repulsive barrier at distances ∼4 times the particle diameter, which is not predicted by either of the two models (Fig. 4a).Proteins with Y p = 50 MPa 100 MPa produce potential wells of similar depth, with Y p = 50 MPa having smaller repulsive barriers (Fig. 4b).Proteins with Y p = 200 MPa, however, show large repulsive barriers and shallow potential wells (Figs.4a and 4b).The difference can be explained by considering the interplay between the energetic and entropic effects.Stiffer proteins induce larger curvatures (Fig. 2), which in turn result in larger curvature-mediated repulsions (Eq.1).On the other hand, the attractive fluctuation-mediated interactions only arise when the protein is stiff enough to locally stifle membrane fluctuations.Thus, there should be a window of optimal protein stiffness.Fig. 4b supports this argument, suggesting proteins with Y p = 100 MPa to be closest to optimal among the cases.

Dynamics of protein aggregation
In all cases investigated here, attractive membranemediated interactions are rather weak compared to the thermal energy, kT (Fig. 4), consistent with the prediction of the continuum model [7].In order to drive the formation of stable protein clusters, cooperative, and pos-sibly multi-body effects are necessary.We investigating these effects by, examining the large-scale dynamics of protein clusters.Among the significant body of work on thermodynamics and statistical mechanics of particle clusters and the aggregation kinetics (see for example [54] and the references within), we have chosen the kinetic theory starting from the following reaction [55][56][57], where C k denotes a cluster of k particles, and K ij and F ij are the corresponding coagulation (formation) and fragmentation (break-up) rates.Assuming detailed balance, as well as a combinatorial model of the number of ways in which particles can form a cluster of a given size, this model leads to an equilibrium cluster size distribution, while allowing for the formation of one macroscopic cluster (details in Methods section "Cluster size distribution").We analyze 3 ms of simulation data, presenting an abundance of these formation/break-up events (Fig. 5a), and obtain dynamical cluster size information from trajectories (details in Methods section "Cluster size distribution").
The most probable cluster size distribution predicted by the kinetic theory (Eq.( 14)) fits very well to our simulation results (Fig. 5d).The µ ′ values serve as a dimensionless chemical potential for exchange of particles between clusters, with lower values hinting at more sta-  ble aggregates.We observe that the increase in protein stiffness disrupts the formation of larger clusters.The maximum cluster size is also achieved for the most flexible proteins (Fig. 5c).
The macroscopic cluster predicted by the theoretical model is highlighted by the g * parameter in Eq. ( 14) denoting the fraction of particles comprising one large cluster.Values directly obtained from the simulation (g * sim ), and values resulting from fits of Eq. 14 to the cluster size distribution (Fig. 5d), show good agreements for Y p = 50 MPa and 200 MPa, while for Y p = 100 MPa, there is an apparent mismatch (ratios cited in Fig. 5d).The probable cause for this is that the cluster growth follows a rather slow dynamics, rendering the formation of macroscopic clusters rare events.The system with Y p = 50 MPa proteins shows the fastest dynamics (compare equilibration timescales in Fig. 5b) explaining convenient formation of a large cluster (Fig. 5c and the ratio g * /g * sim ≈ 1).On the other hand, Y p = 100 MPa proteins have the slowest aggregation dynamics (Fig. 5b), suggesting that the simulation time has not been long enough for these proteins to form the predicted macroscopic cluster.
For Y p = 100 MPa, we have obtained a value of g * = 0.06, meaning that potentially a large cluster incorporating about 6 % of all particles is possible in the simulated conditions.In general, the possibility of the formation of this large cluster signals a phase transition similar to the sol-gel transition in polymers [55].

Free energy landscape of protein aggregation
Free energy landscape of the aggregation process offers the best outlook on the stability of protein clusters and the suggested phase transition.To this end, we have performed equilibrium simulations with different surface concentration of peripheral proteins.We use a reaction coordinate q, corresponding to the continuous dispersion ↔ aggregation transition.We investigate the free energy landscape along this reaction coordinate using histogram reweighting (details in Methods section "Free energy calculation").
The free energy landscape strongly depends on the protein stiffness (Figs.6a and 6b).Increasing the surface concentration results in the free energy minima generally shifting towards more dispersed states (Fig. 6a).But this shift is very small for the Y p = 50 MPa protein, and more drastic for the Y p = 200 MPa.
At the low surface concentration of Γ 1 = 1.5 nmol m −2 , free energy landscape of proteins with Y p = 50 MPa have a clear double-well shape (Figs.6b) pointing to a phase co-existence between the intermediate state "1" and the more aggregated state "2".The same could hold for Y p = 100 MPa proteins, but the results in this case suffer from lack of sampling of the more aggregated state, similar to the observation in the previous section.
At the higher surface concentration of Γ 2 = 9.5 nmol m −2 proteins with the stiffness of 50 and 100 MPa have kept an intermediate state "3", while 200 MPa proteins can only form the more dispersed phase "4" (Figs.6b and 6c).
Free energy landscape along q and protein stiffness clearly shows how increasing the protein stiffness results in a transition towards more dispersed phases (Fig. 6d).We are thus able to verify our hypothesis on the existence of an optimal protein stiffness which minimizes free energy for a given degree of aggregation (Fig. 6d).

Stability of large protein clusters
We have obtained the free energy landscapes in Fig. 6 using equilibrium simulations, which have the downside of limited sampling if large energy barriers exists.In order to investigate if macroscopic clusters of proteins could be stable for long times, we have performed a simulation with an already formed cluster of Y = 200 MPa proteins on the surface of the membrane (Fig. 7a).Even after 1.35 ms, the macroscopic cluster is not dispersed (Fig. 7a).While according to our results on aggregation dynamics, proteins with Y = 200 MPa have the least tendency to form macroscopic clusters (Fig. 5c and 5d).Also, we measured the smallest attractive membranemediated interactions between these proteins (Fig. 4a).
This apparently contradictory behavior can be explained by looking at the simulation snapshots (Fig. 7a and Supplementary Movie 2).The initial protein cluster rapidly induces an endocytic pit.Proteins residing in this pit experience a stabilizing effect due to the fact that escaping it requires passing an unfavorable region of opposing curvature.
To better highlight this stabilizing effect, we compare the root mean squared distance (RMSD) of protein particles from the center of the membrane with membrane particles of a simulation without an initially large protein cluster (Fig. 7a).While for dispersed membrane particles the expected diffusion results in RMSD increasing with time, the protein cluster tends to preserve the compact configuration (Fig. 7a).

Membrane remodeling
Finally, using steered simulation, we look at the free energy landscape of the process in which peripheral proteins aggregate, form an endocytic pit, and finally invaginate the membrane (Fig. 7b).We calculate free energies from biased trajectories via thermodynamic integration (details in Methods section "Free energy calculation").
The free energy landscape reveals a very large barrier facing the pit formation, which explains why this transition is not observed in equilibrium simulations.But once the pit is formed, it is highly stable (Fig. 7b), as we 9.75 10.00   15)) and surface concentration, Γ.(b) Top: Free energies on cross-sections Γ 1 and Γ 2 shown in (a) (gray dashed lines).Shaded regions coincide with numbered states in (a).Bottom: correlation between the inferred bias, ξ, and surface concentration of proteins (Methods section "Free energy calculation") (c) The reaction coordinate q * that minimizes the free energy at each value of Γ.(d) Top left: free energy landscape along q and protein stiffness, Y p .Top right: correlation between the inferred bias and the protein stiffness.Bottom: Free energies on the cross-sections taken at q 1 and q 2 , corresponding to dispersed and aggregated states (color-coded dashed lines in the top plot).showing the evolution of a circular cluster of peripheral proteins.Membrane particles have been rendered transparent to improve visibility.Red particles represent proteins with the stiffness of 200 MPa (Supplementary Movie 2).Bottom: evolution of the root mean squared radial distance (RMSD) of particles, measured from the center of the membrane patch, and normalized with respect to their respective initial value.Results are given for either membrane or protein-supporting particles, where both were initially selected from a similar central disk.(b) Free energy landscape of membrane tubulation resulting from cooperative action of membrane-bending peripheral proteins.Snapshots are from steered simulations in which proteins are forced towards the center of patch (Supplementary Movie 3).also observed in the equilibrium simulation starting with a macroscopic cluster (Fig. 7a).This very stable phase corresponds to a hemispherical pit lined with proteins of a matching curvature.Growth of the pit into a deeper invagination requires incorporation of proteins into a cylindrical wall, which explains the barrier between the pit and the invaginated state (Fig. 7b).Once this barrier is overcome, further growth of the tubule is possible with less energy expenditure.

DISCUSSION
We have presented a dynamical model of biomembranes with flexible membrane-bound proteins to study their aggregation behavior.The model reproduces intended membrane curvatures, as well as realistic membrane viscosity and in-plane diffusion of proteins.It is thus suitable for studying thermodynamics and kinetics of membranes and membrane-associated proteins, especially in large-scale simulations elucidating collaborative actions and macroscopic phenomena.
Our model is purposefully void of a direct interaction between protein pairs, leading us to conclude that the measured pairwise interactions reveal membranemediated attractive forces.Higher protein stiffness, which in our model also includes stronger membranebinding, results in larger fluctuation-based attraction, but also enhances repulsions originating from membrane bending energy.This contrariety results a selective mechanism for protein flexibility.This selectiveness was further highlighted when we discovered that more flexible proteins are better at forming large clusters.
We obtained an optimal stiffness value for a given aggregation tendency (Fig. 6d).The local membrane curvature in our model that would result from proteins with this optimal stiffness of 63 MPa, is 0.032 nm −1 , which is in very good agreement with the molecular dynamics simulation of STxB yielding curvatures 0.034 to 0.035 nm −1 [41].
We generally observe the magnitude of attractive membrane-mediated interactions to be small.Yet the collaborative effect of these interactions seem to suffice to form loosely-bound clusters (Fig. 5a).These clusters can still impact the macroscopic membrane conformation when they form pits and experience a subsequent stabilizing effect (Fig. 7).We uncovered a rather complex pathway and large energy barriers in-volved in formation and stabilization of macroscopic protein clusters and membrane deformations.We conclude that the suggested pathway of spontaneous aggregation leading to membrane tubulation may need more elaboration.The high energy barriers suggest that proteins are more probably attracted to the already existing membrane conformations.
The rare events indicated in membrane remodeling pathway deserve more detailed analysis in the future.The approach laid out here paves the way for direct comparison of simulation and experiment to study processes such as clathrin-independent endocytosis of toxins or signaling molecules.It also offers a quantitative tool for designing finely controlled drug delivery vehicles that rely on membrane mechanics for cellular entry.

Elastic model of membrane and peripheral proteins
For a homogeneous membrane the elastic energy is very well described by the Helfrich model [58][59][60][61], in which H (x) and G (x) are the mean and the Gaussian curvatures of the membrane at position x, and κ and κ are the corresponding bending rigidity and Gaussian curvature moduli.H 0 (x) represents a spontaneous curvature, arising from asymmetries between the two leaflets and the shapes of individual lipid molecules.
To build an elastic model of the peripheral protein, we assume a "preferred" intrinsic curvature, c 0 , for the protein, such that the bound protein is unstrained when the underlying membrane patch acquires the mean inward curvature of c 0 .But the flexibility of the protein as well as the protein-membrane coupling means that in mechanical equilibrium, the protein has non-zero principal curvatures c ′ 1 and c ′ 2 .The curvatures of the protein and the membrane are complementary in the sense that c ′ 1 + c 1 = c 0 and c ′ 2 + c 2 = c 0 , with c 1 and c 2 being principal curvatures of the membrane (Supplementary Fig. 1).We employ the Kirchhof-Love small-deformation plate theory as a continuum model of the protein, resulting in the protein strain energy per unit area, in which Y p is the effective Young modulus of the protein, ν p is its Poisson's ratio and d p is its thickness (Supplementary Fig. 1).Having the energy densities of both the protein and the membrane, the total energy density for which the particle-based model will be parametrized is, where Ap Am denotes the surface area ratio between the peripheral protein and the membrane in the binding region, and ✶ p (x) is an indicator function equal to one where a peripheral protein is bound to the membrane, and zero otherwise.

Simulations
Details and parameters of the force field for the unmasked and masked regions are given in Supplementary Information.We have used over-damped Langevin dynamics with anisotropic diffusion tensors to describe the motion of membrane and protein particles [38,39].We have considered hydrodynamic coupling to an aqueous solvent with the viscosity of η = 0.890 mPa s, and have included hydrodynamic interactions between nearest-neighbor particles (the "Hydro NN" model described in ref. [39]).All simulations have been carried out at the biological temperature of T = 310 K.We have coupled the in-plane degrees of freedom to the Langevin piston barostat to simulate tensionless membranes [62].
The in-plane fluidity of the membrane is modeled using bond-flipping Monte Carlo moves [37].One move constitutes switching a bond shared between two adjacent bonded triangles to the crossing diagonal.We evaluate the acceptance probability of a proposed bondflip or force field switch, which is also implemented via Monte Carlo moves, via the Metropolis-Hastings algorithm [46].In order to guarantee that the masking of the forcefield does not disturb bond-flipping, we have separated the masking of the force field into two distinct contributions, one that preserves the stiffness of in-plane bonds, and an additional stiffening potential (details in Supplementary Information).Bond-flipping is decided based on the former set of interactions, and results in the force field parameters of flipped bonds to reset to the original unmasked values.We have monitored the potential energy of bonded interactions throughout the simulations to ensure that this process does not result in any unwanted energy fluctuations.
To measure the effective local curvature around each protein-supporting particle (Fig. 2), we have fitted a bivariate quadratic function to the position of each particle and its nearest neighbors.These spontaneous values of curvatures have been sampled for a random selection of particles for each frame of the trajectory.

Statistical mechanics of liquid mixtures
We assume rotational isotropy in the 2D fluid, and consider the pair-correlation function between the species α and β as g αβ (r).Defining h αβ (r) = g αβ (r) − 1, the so-called direct correlation function, c αβ (r), is defined implicitly [63], where Γ's denote surface densities.Eq. ( 7) translates to where H and C are 2 × 2 matrices with elements Hαβ = (Γ α Γ β ) Eq. ( 8) is formally known as the Ornstein-Zernike relation [63,64].Several "closure" relations are available which tie h αβ , c αβ , and g αβ via a pairwise potential u αβ between α and β.We have employed the convoluted hypernetted-chain (CHNC) relation [65], The procedure for obtaining the potentials u mm , u mp , and most importantly, u pp , is as follows: we calculate discretized estimates of the three pair correlation functions, g mm (r), g mp (r), and g pp (r), from the simulation trajectories, and substitute numerically obtained 2D Fourier transforms into Eq.8 to obtain c αβ 's and subsequently u αβ 's from Eq. ( 9).

Cluster size distribution
A necessary condition for the existence of an equilibrium cluster size distribution is the following detailed balance condition [55,56], where at this point, Λ and a k 's are arbitrary set of parameters, for which different interpretations exist.One popular approach is the combinatorial interpretation of k! a k as the number of ways that a cluster of size k can be formed from k distinguishable particles [55][56][57].If n = (n 1 , n 2 , . . ., n M ) is considered as the state vector, with n k representing number of clusters of size k, the equilibrium distribution is given as [55], where A is the projected surface area, and Z (Λ, A, M ) is the normalizing factor i.e. the clustering partition function.The most probable cluster size distribution is thus achieved when, in which the presence of one large cluster of size k 0 is allowed.If this largest cluster incorporates a fraction g * of particles, the mass conservation condition reads as, Based on the aforementioned combinatorial interpretations, we have assumed a k ∝ c k , with c being a constant.The justification is that the number of ways that a new particle can be added to a 2D cluster of proteins depends on a maximum coordination number.This assumption leads to the following expression for the most probable cluster size distribution, where µ ′ = µ − ln c has been substituted.To obtain cluster size information from trajectories, we have used the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm [66].This method can discover point clusters based on a pairwise distance metric and a tolerance measure for particles belonging to a cluster.We have set the tolerance value to be twice the lattice parameter of the model.
The fraction of particles comprising one macroscopic cluster, g * sim , is calculated from the mean of maximum cluster sizes shown in Fig. 5c.

Free energy calculation
We have used the mean of inverse pairwise distances between peripheral proteins as the clustering reaction coordinate, In a single simulation, it is trivial to calculate the histogram of microstates along the reaction coordinate q.
To correctly combine results from different simulations, we have employed the multiple histogram method developed by Ferrenberg and Swenden, which is the precursor to the well-known Weighted Histogram Analysis Method (WHAM) [67,68].We assume that each simulation, in which a different number of particles have been tagged as peripheral proteins, is a biased simulation, corresponding to a Hamiltonian of the form H = H 0 +ξ q.If this assumption holds, the bias parameter ξ should correlate in the two studied cases with (i) the surface concentration, Γ, or (ii) stiffness of proteins, Y p .We have estimated the value of ξ for each simulation via comparing H between all simulations.We observe that the expected correlations exists between ξ and Γ or ξ and Y p (Figs. 6b and 6d), ratifying our approach for inferring the bias from simulations.
To obtain the free energy landscape of membrane tubulation process (Fig. 7b), we have set up 10 simulations, in which particles supporting peripheral proteins are pulled toward the center of the box with a fixed force.Free energy difference is calculated using the thermodynamic integration, where we have used the sum of radial coordinates of forced particle as the parameter λ.

Software
Simulations based on the particle-based membrane model [37][38][39] are performed using an in-house specific-purpose C++ software.Python libraries Numpy, SciPy and Matplotlib are used for analysis and presentation of the results [69][70][71].The software package Visual Molecular Dynamics (VMD) is used for some visualizations [72].

SOFTWARE AVAILABILITY
All the in-house developed software used in this study are available from the corresponding authors upon request.

− 1 ]c 1 YFigure 2 :
Figure 2: Local membrane curvature in the vicinity of peripheral proteins with the given stiffnesses.Distribution of instantaneous curvature values is shown on the right.Dashed lines pinpoint the mean of each distribution, while the shaded regions around the time series correspond to one standard deviation.The arrows point to the estimated curvatures from the parameter-space optimization procedure, obtained from energy distributions in curvature space (Fig.1).

Figure 3 :
Figure 3: Kinetics of membrane and peripheral proteins: (a) Mean-squared displacement for in-plane motion of proteins.Dashed lines represent power-law fits to the shown range of lag times.Values of the corresponding exponents of time t and the effective long-time diffusion coefficients are color-coded for the given protein stiffnesses.(b) Time-correlation functions of the in-plane shear stress of the membrane.Values of membrane surface viscosity, η m are obtained using the Green-Kubo relation and are given in each case in units of surface Poise (sp).

Figure 4 :
Figure 4: The membrane-mediated interaction between peripheral proteins: (a) Interaction potential as a function of inter-particle distance, for proteins with the given stiffnesses.Results are given for simulations with the protein surface concentration of Γ = 2.0 nmol m −2 .Predictions of the theoretical model for pairwise curvature-mediated (repulsive) interactions, U rep , as well as the fluctuation-mediated (attractive) interactions, U att , are plotted for comparison (Eqs.(1) and (2)).Contact angles (θ p ) are calculated based on the measured mean curvatures given in Fig.2.(b) Effect of the protein stiffness on the shape and strength of the interaction potential.Well-depth, ∆u w , and barrier height, ∆u b , are measured according to the inset schematic shown in (a).The red and gray shaded regions in (a) are where the local minima and maxima of potential energies have been sampled.

Figure 5 :
Figure 5: Dynamics of protein clusters on the surface of the membrane: (a) Simulation snapshots showing the formation and break-up of peripheral protein clusters for the protein surface concentration of Γ = 2.0 nmol m −2 (Supplementary Movie 1) (b) Time evolution of the weighted mean cluster size (number of particles in a cluster) for peripheral proteins with the given stiffnesses.Dashed lines are fits of an exponential saturation model N t = N ∞ + ( N 0 − N ∞ ) exp (−t/τ ), with values of timescale τ , and equilibrium mean cluster size N ∞ , given in each case.(c) Maximum cluster size achieved during the simulation.The background of each plot shows the evolution of the histogram of cluster sizes.(d) Equilibrium distribution of cluster size.Dashed lines are fits of Eq. (14) with the best fit parameters and the R 2 score of the fit given on the plots.

Figure 6 :
Figure6: Free energy landscape of protein aggregation: (a) Free energy landscape along the aggregation reaction coordinate, q (Eq.(15)) and surface concentration, Γ.(b) Top: Free energies on cross-sections Γ 1 and Γ 2 shown in (a) (gray dashed lines).Shaded regions coincide with numbered states in (a).Bottom: correlation between the inferred bias, ξ, and surface concentration of proteins (Methods section "Free energy calculation") (c) The reaction coordinate q * that minimizes the free energy at each value of Γ.(d) Top left: free energy landscape along q and protein stiffness, Y p .Top right: correlation between the inferred bias and the protein stiffness.Bottom: Free energies on the cross-sections taken at q 1 and q 2 , corresponding to dispersed and aggregated states (color-coded dashed lines in the top plot).

Figure 7 :
Figure 7: Stability of large protein clusters and the resulting membrane remodeling: (a) Top: simulation snapshotsshowing the evolution of a circular cluster of peripheral proteins.Membrane particles have been rendered transparent to improve visibility.Red particles represent proteins with the stiffness of 200 MPa (Supplementary Movie 2).Bottom: evolution of the root mean squared radial distance (RMSD) of particles, measured from the center of the membrane patch, and normalized with respect to their respective initial value.Results are given for either membrane or protein-supporting particles, where both were initially selected from a similar central disk.(b) Free energy landscape of membrane tubulation resulting from cooperative action of membrane-bending peripheral proteins.Snapshots are from steered simulations in which proteins are forced towards the center of patch (Supplementary Movie 3).