Protein crowding and lipid complexity influence the nanoscale dynamic organization of ion channels in cell membranes

Cell membranes are crowded and complex environments. To investigate the effect of protein-lipid interactions on dynamic organization in mammalian cell membranes, we have performed coarse-grained molecular dynamics simulations containing >100 copies of an inwardly rectifying potassium (Kir) channel which forms specific interactions with the regulatory lipid phosphatidylinositol 4,5-bisphosphate (PIP2). The tendency of protein molecules to cluster has the effect of organizing the membrane into dynamic compartments. At the same time, the diversity of lipids present has a marked effect on the clustering behavior of ion channels. Sub-diffusion of proteins and lipids is observed. Protein crowding alters the sub-diffusive behavior of proteins and lipids such as PIP2 which interact tightly with Kir channels. Protein crowding also affects bilayer properties, such as membrane undulations and bending rigidity, in a PIP2-dependent manner. This interplay between the diffusion and the dynamic organization of Kir channels may have important implications for channel function.

Cell membranes are crowded and complex environments. To investigate the effect of protein-lipid interactions on dynamic organization in mammalian cell membranes, we have performed coarsegrained molecular dynamics simulations containing >100 copies of an inwardly rectifying potassium (Kir) channel which forms specific interactions with the regulatory lipid phosphatidylinositol 4,5-bisphosphate (PIP 2 ). The tendency of protein molecules to cluster has the effect of organizing the membrane into dynamic compartments. At the same time, the diversity of lipids present has a marked effect on the clustering behavior of ion channels. Sub-diffusion of proteins and lipids is observed. Protein crowding alters the sub-diffusive behavior of proteins and lipids such as PIP 2 which interact tightly with Kir channels. Protein crowding also affects bilayer properties, such as membrane undulations and bending rigidity, in a PIP 2 -dependent manner. This interplay between the diffusion and the dynamic organization of Kir channels may have important implications for channel function.
Cell membranes are crowded and complex environments, containing up to 50% protein by mass, which is equivalent to ca. 25% protein by cross-sectional area 1 . They contain a multiplicity of lipid species distributed asymmetrically between the outer and inner leaflets of the bilayer [2][3][4] . Due to the high area fraction of protein in cell membranes, protein crowding occurs. This has been implicated in a number of biological processes [5][6][7] , and protein nanoclustering has been shown to be an important factor in membrane organisation 8,9 . It has been shown that membrane proteins co-diffuse as complexes with lipids 10 , and that protein crowding affects both protein and lipid diffusion 11-13 . There is a growing appreciation of the importance of specific protein-lipid interactions and their role in both regulating protein activity and contributing to the dynamic high order organization of cell membranes [14][15][16] . Molecular simulations in particular have proved a useful tool to study lipid interactions with membrane proteins [17][18][19] .
In addition to probing the nature of protein-lipid interactions, there has been much interest in how complex lipid mixtures may generate dynamic organization in lipid bilayers. This has embraced e.g. lateral nanodomains of specific lipid composition 20,21 , and also the contribution of lipid properties to membrane curvature 22 and fluctuations 23 . Despite impressive progress 7 much remains to be understood regarding how the molecular interactions between protein and lipids, including crowding, influence higher-level membrane organization and dynamics 24 . Advances in both coarse-grained (CG) force fields 25 and in computational power mean that large-scale CG molecular dynamics (MD) simulations are increasingly being used to provide insights that can be directly linked to experiments 26,27 .
In order to investigate the role of protein crowding and lipid complexity on the organization and dynamics of a model membrane containing a protein whose activity is regulated by interactions with lipids, we have focused on the ion channel Kir2.2. This is a member of the inward rectifier family of potassium channels, which have well-characterized and specific lipid interactions [28][29][30] . Kir channels help to restore the electrical voltage across a membrane to its resting potential in excitable cells 31 . Kir2.2 channels are tetrameric proteins containing a classical potassium channel transmembrane domain coupled to a large cytoplasmic domain 32 . Kir channels are activated by phosphatidylinositol 4,5-bisphosphate (PIP 2 ), an anionic lipid present in the inner leaflet of mammalian plasma membranes. PIP 2 binds to four identical interaction sites on the surface of the tetrameric channel protein 32,33 . Secondary anionic lipid species (e.g. phosphatidyl serine) may further modulate the activity of the channel, and an additional anionic interaction site on each subunit has been identified 29,34,35 . Additionally, cholesterol inhibits the channel via direct protein-lipid interactions 36,37 . The cytoplasmic domain of the homologous Kir2.1 channel has been shown to form clusters in the presence of channel-organising protein, PSD-95 (see below) 38 .
There is mounting biophysical evidence for ion channel clustering in cell membranes, and for its functional importance e.g. 39 . Clustering of the canonical bacterial potassium channel KcsA is influenced by anionic lipids 40 , with the channel forming large salt-induced complexes 41 . In reconstituted membranes KcsA can form extended clusters up to a micrometer in size which are observable by confocal microscopy 42 . Other bacterial channels have also been demonstrated to cluster, for example including the mechanosensitive channel MscL 43 . Several studies have implicated lipids in channel clustering mechanisms, e.g. lipid microdomains were identified as being important in Ca v 2.1 clustering 44 . Early functional studies of A-type potassium channels suggested the formation of clusters of ~50 channels 45 , and fluorescence imaging studies revealed that L-type Ca 2+ channels form clusters of ~100 nm diameter 46 . The canonical MAGuK (membrane-associated guanylate kinase) scaffold protein PSD-95 aids the formation of potassium channel clusters on postsynaptic membranes 47,48 . A number of studies have indicated that Kir channels interact with PSD-95 [49][50][51][52][53] and the isolated cytoplasmic domain of Kir2.1 has been shown to form clusters on association with PSD-95 38 . It is therefore timely to examine the extent to which Kir channels cluster when embedded in model membranes at a surface density comparable to that in mammalian cell membranes 1 , and how the resultant dynamic organization of channel molecules may be influenced by the lipid environment.
Here we present large-scale CG molecular dynamics simulations of Kir2.2 channels embedded in a range of membranes of varying degrees of lipid complexity and protein crowding. We dissect the interplay of protein crowding and lipid dynamics, and show that the interactions between Kir channels and PIP 2 play a key role in the dynamic organization of a biologically realistic model of cell membranes.

Results and Discussion
Simulations of Crowded Kir-Containing Membranes. The simulations performed are summarized in Table 1 (see also Supplementary Table S1). To explore the dynamic organization of Kir channels we simulated 144 copies of the Kir2.2 channel embedded in a plasma membrane (PM) model. This provides a membrane area (of 137 × 137 nm 2 ) with a sufficient density of proteins for clustering to be observed. The model was simulated for 50 μs to allow sufficient time for the protein and lipid dynamics to be extensively explored (Fig. 1a). The bilayer of the PM model contained the major lipid species observed in a mammalian plasma membrane 54 . This model has previously been used to explore GPCR clustering 55 . In order to explore the effect of PIP 2 -Kir channel interactions on the behavior of the ion channel, a 20 μs simulation of 144 Kir2.2 channels in a control plasma membrane lacking PIP 2 was performed (noPIP2; Table 1). Finally, to explore the effects of lipid complexity, an additional control simulation of 144 Kir2.2 channels in a simple PC membrane was performed for 20 μs.
A summary of lipid headgroup interactions with Kir channels observed in the simulations is given in Fig. 1b. In the simulation containing PIP 2 , headgroup interactions with Kir2.2 channels extended further from the transmembrane domain onto the cytoplasmic domain of the channel (comparing top and middle panels in Fig. 1b). The radial distribution functions (Fig. 1c) indicate that in the PM simulation, the principal lipids surrounding Kir2.2 channels are PIP 2 in the cytoplasmic leaflet and GM3 in the outer leaflet. PS was the second most populous inner leaflet lipid in the annular shell around Kir2.2 channels. In noPIP2 simulations, the abundance of PS at the protein surface was increased, but not to the extent for PIP 2 in the PM simulations. Thus it can be seen that the major interactions of the protein are with PIP 2 in the inner leaflet, and with GM3 in the outer leaflet, as has been shown to be the case for the EGF receptor 56 .
Lipid Complexity and Dynamic Channel Protein Organization. The organization of Kir2.2 channels after 20-50 μs of simulation in membranes of differing lipid complexity are illustrated in Fig. 1a. Substantial differences in the clustering of Kir2.2 channels are apparent. In membranes with complex lipid mixtures (PM and noPIP2), the ion channel clusters are smaller and appear to be more branched (especially in the PM simulation). If we track cluster formation (Fig. 2) and cluster size (Fig. 3a) as functions of time this difference is especially evident. After 15 µs most of the channel proteins in the PC membrane form a single cluster containing ca. 120 proteins ( Fig. 2c) whereas at the same time in the PM simulation clusters range up to ca. 30 proteins in size (Fig. 2a), with the noPIP2 simulation exhibiting intermediate behavior (Fig. 2b). It is evident from the average cluster size vs. time that Kir clusters form more rapidly in the PC membrane (Fig. 3a). There is also a noticeable difference in cluster shape between the different bilayer compositions. Clusters in the more complex lipid mixtures are less linear in shape than in the PC membrane. Cluster shape can be characterized by the number of neighbors each protein has: clusters with a linear shape typically contain proteins with only two neighbors, whereas clusters that are more branched contain a high fraction of proteins with more than two neighbors ( Fig. 3b and c). Cluster shapes can be compared by looking at the largest cluster formed of ~100 channels, occurring after 40 μs in PM and 20 μs in the PC simulations. The proportion of proteins with  Table 1 for details). The membrane is viewed from the inner (cytoplasmic) side, with proteins colored according to the cluster size, and lipids colored dark grey. The approximate simulation box size (shown in grey) is 137 × 137 nm 2 . (b) The Kir2.2. structure with lipid headgroup interactions mapped onto the channel surface using a greyscale from white (no interactions) to grey (many interactions). Residues side chains are shown for those residues at which lipid headgroup interactions account for more than 2.5% of the total interactions each. The approximate location of the lipid bilayer is indicated by two broken lines on each structure, with the transmembrane (TM) and intracellular (IC) domains indicated on the diagram corresponding to the PC simulation. (c) Radial distribution functions of the distribution of lipid headgroups around Kir2.2 channels for the PM and noPIP2 simulations.
only one neighbor is similar in both PM and PC simulations, over the timescale that the largest cluster is formed. However, the proportion of proteins with two or more than two neighbors appear to diverge: in the PM simulations, 35-40% of proteins have three or more neighbors, whereas in the PC simulations only ~25% have three or more neighbors.
What underlies this difference in organization? The frequency and pattern of protein-protein interactions, as estimated via the spatial frequency of pairwise interactions (Fig. 3c) and the map of interacting residues (Fig. 4), are similar in all three simulations and are largely mediated by the intracellular domain of Kir. Furthermore, even though Kir2.2 binds PIP 2 molecules close to its intracellular domain (see Fig. 1c), the presence or absence of PIP 2 does not seem to influence channel clustering. Therefore, the influence of lipids on the dynamic organization and clustering of Kir channels does not appear to be mediated directly through protein-lipid interactions, but more likely is an effect of changes in lipid diffusion within the plasma membrane upon increasing lipid complexity.

Lipid Diffusion
Rates. The diffusion rates of the component lipids and proteins in the different membrane simulations appear to correlate with the clustering behavior of the channel proteins (Fig. 5a). Both the channel proteins and all lipid species in PM and noPIP2 simulations diffuse more slowly than in the PC simulation. For example, for the channel proteins the diffusion coefficients are ~2 vs. 4 μs 2 /s α in PM/noPIP2 vs. PC respectively,  To investigate which lipid species is responsible for the reduction in diffusion, further simulations were performed for a variety of complex lipid membranes. These were scaled down in size to yield a system with 9 Kir2.2 channels embedded in a smaller membrane, whilst maintaining the same protein:lipid ratio as for the larger systems (Tables 1 and S1). In order to understand the effect of specific protein-lipid interactions on diffusion, simulations of this smaller system were run for a plasma membrane (PM) model, and for a PM with PIP 2 removed, or GM3 removed, or with both of these lipid species removed. To assess the contribution of cholesterol (Chol) to membrane dynamics, simulations of 9 Kir2.2 channels in a PM with cholesterol removed, or with PIP 2 , GM3 and Chol all removed were also run (see Table 1 for details).
The effect of removing the lipid species which interact with the protein (that is, PIP 2 and GM3 -see above) is to increase the diffusion coefficient slightly for all lipid species (Fig. 5b and Table S2). Removal of either PIP 2 or GM3 mainly alters the diffusion coefficients of the lipids in the same leaflet, (i.e. in the inner or outer leaflet respectively). It is notable that this trend does not hold for PS; even though PS and PIP 2 are both in the inner leaflet of the plasma membrane, on removal of PIP 2 , the diffusion coefficient of PS does not increase as much as that of e.g. PE. This may be because PS retains a high level of interaction with Kir channels. Cholesterol is present at equal levels in both leaflets. Its diffusion coefficient was affected slightly more by the removal of GM3 than that of PIP 2 . Removal of both PIP 2 and GM3 caused a slightly greater increase in the diffusion coefficients of all remaining lipid species than when PIP 2 or GM3 were removed individually, and the effect appeared to be additive, for example POPE diffusion was 10 μs 2 /s α in the PM_s simulation, increasing to 12 μs 2 /s α on removal of PIP 2 (noPIP2_s simulation), and to 11 μs 2 /s α on removal of GM3 (NoG_s simulation). After removal of GM3 and PIP 2 (NoPIP2G_s simulation), POPE had a diffusion coefficient of 13 μs 2 /s α . Similar trends were observed for sphingomyelin, cholesterol and POPC.
The effect of cholesterol on membrane diffusion has been much investigated and has been shown to slow diffusion (see e.g. [57][58][59]. As a control, we removed cholesterol, both as the sole species taken out of the plasma membrane and concurrently with the removal of PIP 2 and GM3. Whereas removal of PIP 2 and/or GM3 caused little increase in the diffusion coefficient of the Kir channels and the remaining lipid species (the maximum increase on removal of PIP 2 and GM3 was for sphingomyelin (Sph) which increased by 5 μs 2 /s α from PM_s to noP2G_s) in the membrane, removal of cholesterol caused a much greater increase in the diffusion coefficient of all remaining lipid species. Removal of cholesterol from the plasma membrane has a particularly marked effect on sphingomyelin diffusion, as might be anticipated 60 (although the nature of the cholesterol-sphingolipid interaction in living cell membranes is currently the topic of much debate 61 ). Moreover, while the removal of PIP 2 and/or GM3 had little effect on the diffusion coefficient of Kir2.2 channels (an increase from 0.4 μs 2 /s α in PM_s to 0.6 μs 2 /s α in noP2G_s), removal of cholesterol almost tripled the protein diffusion coefficient (from 0.4 μs 2 /s α in PM_s to 1.1 μs 2 /s α in noChol_s). However, even with the removal of much of the lipid complexity, protein and lipid diffusion is not increased to the levels observed in the single-component PC membrane. Thus, diffusion of the lipid PC is 3 μs 2 /s α lower in PM simulations (comprised of a lipid mixture containing at least PC, PE, PS and sphingomyelin) than in PC simulations (which contain only PC).

Protein Crowding Decreases Diffusion Coefficient.
In order to investigate the effect of protein crowding on diffusion and channel clustering, a simulation was performed in which only 36 Kir channels were embedded in a 120 × 120 nm 2 membrane, thus giving a system (PMsparse; Table 1) with 25% of the protein density of the PM simulations (Fig. 6). Protein clusters in this sparsely populated membrane were much smaller, with typically at most dimers of Kir channels formed on the timescales simulated here. As a further control, a simulation of the plasma membrane containing no protein (PMnoprot) was also performed. Decreasing protein crowding caused both protein and lipid diffusion coefficients to increase significantly, although not quite to the values for the PMnoprot system. Thus, reducing the protein density to the extent that crowding appeared to no longer play a role (PMsparse) increased lipid diffusion coefficients to values of ~60% of the PMnoprot coefficients (Fig. 6c). A further control simulation of a PC membrane containing no protein (PCnoprot) also showed lipid diffusion to increase greatly on removal of protein ( Figure S1A).

Anomalous Diffusion Behavior is Modulated by Crowding. Protein crowding had a marked effect
on the diffusive behavior of the protein (Fig. 7). In all of the crowded simulations (i.e. PM, noPIP2, PC) the value of the anomalous exponent α for proteins initially decreases from 0.7-0.8 at Δt of 1-10 ns to a minimum value of 0.5-0.6 at Δt in the range 10-100 ns. The value of α then increases sharply when measured over timescales of hundreds of nanoseconds, reaching ~1.0 at a timescale of 10 μs. As anticipated, proteins in the more complex lipid mixture simulations (PM and noPIP2) have overall lower anomalous exponents (and are therefore more sub-diffusive) than those of the simple PC membrane. Similar protein behavior is observed in the small simulations ( Figure S2). At large timescales (Δt > 1 μs) the protein α values of the PC membrane system start to exceed 1, which may be due to the large protein clusters, in which the proteins tend to move in similar direction, and to a lack of sampling at larger timescales.
When crowding is decreased, as in the PMsparse simulation, the magnitude of the anomalous diffusion coefficient for the protein is higher when measured over timescales <1 μs compared to the crowded PM, and the variation in α as the measurement timescale is varied becomes much less pronounced (the protein α values remain between 0.7 and 0.8 for Δt <1 μs). At larger (Δt) timescales (1-10 μs), sub-diffusive behavior of the PMsparse system becomes more like that of the crowded PM system.  Table S2. (b) Diffusion coefficients for the small simulations systems, again with protein diffusion coefficients as an inset (see Table 1 for details of system compositions, and Table S2 for diffusion coefficients listed). The small system size is indicated by the inset, as before with the black scale bar showing 30 nm. (c) Color codes for each simulation in A and B are shown along with lipid proportions for each simulation system, with lipid coloured: PC (black), PE (grey), PS (blue), PIP 2 (orange-yellow), Sph (green), GM3 (magenta) and cholesterol (cyan).
By varying the concentration of the Kir channel in the membrane, the effect of protein-lipid interactions on lipid diffusion can also be inferred (Fig. 8). At high levels of crowding, when almost all the PIP 2 and GM3 lipids are engaged in protein-lipid interactions, PIP 2 and GM3 have much lower α values than other lipid species at low Δt timescales (in the order of 10-100 ns), indeed sub-diffusive behavior is similar to that of the protein (Fig. 8   Figure 6. Snapshots of the (a) uncrowded PMsparse and (b) crowded PM systems, both at 20 μs, with proteins colored according to cluster size (as in Fig. 1a), viewed from the cytoplasmic side. (c) Effect of crowding on diffusion, shown as diffusion coefficients of proteins and each lipid species in the crowded PM (orange-yellow), the uncrowded PMsparse (lime green) and the PMnoprot (red) simulations. Diffusion coefficients are also listed in Table S2. and Figures S1, S2 and S3). However, if protein concentration is decreased (in the PMsparse simulation), the sub-diffusive behavior of PIP 2 becomes similar to that of other lipid species. In simulations without protein, PIP 2 α values are indistinguishable from those of other lipid species. GM3, however, continues to have much lower α values than other lipids when proteins are removed. This may be due to GM3 clusters forming in the membrane without the presence of protein 54 , causing marked sub-diffusive behavior of GM3 at the 10-100 ns timescales.
We have observed a complex interplay between the clustering of the Kir channel protein and its diffusive dynamics. In particular, the clustering of proteins in the membrane, driven by crowding, leads to markedly anomalous (and slowed) dynamics of those lipids, which interact most strongly with the channel protein, namely PIP 2 and GM3. In turn, slower diffusion reduces the rate of protein clustering, resulting in more compact and branched protein clusters.

Membrane Undulations.
In previous simulations of a protein-free PM, the membrane was observed to undergo substantial vertical fluctuations leading to undulatory motions, reflecting a significantly reduced bending rigidity of the membrane relative to a simple PC bilayer 23,54 . Such large movements (ca.±15 nm along the bilayer normal) are also seen in the PMsparse simulation (Fig. 9). In contrast, the fluctuations are smaller (ca.±5 nm along the bilayer normal) in the PC and noPIP2 simulations. Interestingly, in the PM membrane (containing 144 Kir molecules) the magnitude of the fluctuations is intermediate between these two extremes. The overall pattern of the magnitude of bilayer undulations is PMsparse >PM >noPIP2 ≈ PC (Fig. 9b).
Simulations of membranes lacking PIP 2 displayed little change in membrane shape during the course of the simulations (Fig. 9a), with the membrane remaining planar. However, in PM simulations, where compositions were identical except for the presence of PIP 2 , changes in membrane shape were observed, with substantive undulations in the membrane appearing. This propensity for membrane deformation increased as the density of protein, and hence crowding, decreased, indicating that in the PM system a degree of protein crowding and clustering flattened the membrane. The observed effect of PIP 2 conflicts with studies of simple model membranes formed by e.g. amphiphiles 62 which indicate that increasing the surface charge on a membrane lead to an increase in bending rigidity. We presume this difference arises because our model membranes do not form a symmetric bilayer and also because their lipids cluster laterally. Thus, both lipid composition (the presence of anionic lipids) and clustering of membrane proteins can modify the bending rigidity of a model membrane.
Membrane undulations were quantified by calculating the power spectra of the height fluctuations of the bilayer (although the asymmetry of the bilayers precludes calculation of the bending modulus from these Scientific REPORTS | 7: 16647 | DOI:10.1038/s41598-017-16865-6 spectra 23,63 ). The power spectra show that for low wavenumbers (q < 0.1 nm −1 , corresponding to undulations in the range 10-20 nm) there was an effect of lipid composition on the intensity of undulations such that the magnitudes of the undulations followed the order PM >noPIP2 ≈ PC (Fig. 9c). At higher wavenumbers (q > 0.5 nm −1 , or undulations <2 nm) the power spectra converge, indicating little effect of the different lipid compositions. Thus, the effect of adding of PIP 2 affects the dynamics of the membrane on a length scale of 5 nm and above, corresponding approximately to the size of a Kir channel molecule.

Implications.
We have revealed a link between the clustering of Kir channel proteins and the anomalous diffusion of both proteins and lipids. Anomalous diffusion has been observed in a number of previous simulation studies of membrane proteins 11,[64][65][66] , and in particular the protein diffusion coefficients observed here agree well with e.g. membrane protein diffusion rates in PC/PG liposomes measured using fluorescence correlation spectroscopy 67 . Our simulations allowed us to analyze the effect of protein crowding and clustering on lipids that interact specifically with the channel protein (PIP 2 and GM3), and which, in the case of PIP 2 is an allosteric activator of the channel 68 . Furthermore, our simulations provide an example of a 'clustered protein lattice' in a reasonably complex model of a cell membrane, and reveal how such cluster formation may influence both lipid and protein dynamics within a membrane.
Of course, the simulations described remain a simplification of real cell membranes. In particular, the clustering of Kir channels in the simulations is brought about purely by crowding (albeit at 'physiological' densities) and other proteins are not present. It is likely that the cytoskeleton immobilizes certain membrane proteins [69][70][71] , an effect which in the case of Kir channels may be mediated in postsynaptic membranes via PSD-95 38 . Our simulations suggest that such cytoskeletally-mediated clustering may build upon an intrinsic dynamic organization of Kir channels which is brought about by crowding and diffusional restriction in complex cell membranes 72 .
It will be important to understand the impact of channel clustering on Kir channel activity. Given the interactions of (clustered) Kir channels with the activating lipid PIP 2 , it will be of interest to compare our observations in what remains a relatively simple membrane model with behavior in more complex laterally anisotropic cell membranes containing e.g. cholesterol-induced nanodomains. We are also not able to observe large-scale protein conformational changes using CG MD (for instance, those related to channel activation). The gating kinetics of Kir channels, and how this may effect, and be effected by, the crowding of Kir channels is an interesting further question, which will require further work, although a recent communication suggests a link in the case of KcsA 73 . Furthermore, understanding the localization and behavior of complexes of Kir and other channels is directly relevant to a number of disease states, e.g. arrythmogenic cardiac disease 74 . The first concerns simulation size effects: diffusion is expected to increase with system size 75,76 , although this is not seen in the simulations presented here. This is perhaps due to there being multiple copies of Kir channels even in the 'small' systems and the formation of compartments due to protein clustering in both 'small' and 'large' systems.
There also may be issues of convergence and of the behaviour of the coarse-grained forcefield in large scale simulations with many interacting proteins 77,78 , and in the forcefield representation of certain lipids, eg. GM3 79 . However, we should note that contrary to some expectations, the MARTINI forcefield seems if anything to underestimate the strength of protein/protein interactions within membranes once a rigorous comparison with experimental data is undertaken 80 . We note that in this cited study a different way of accounting for electrostatics is employed from in the current study (PME rather than the switch/cutoff function used here). However, we note that a number of other studies using the switch/cutoff representation of electrostatics have used the MARTINI force field to specifically identify experimentally validated protein-protein interfaces [81][82][83][84] . Furthermore, a recent study of the interaction of gangliosides with the WALP transmembrane helix peptide and with aquaporin, showed that, when using polarizable water, PME electrostatics did not have an effect on protein-lipid or lipid-lipid interactions 79 . Whilst the use of polarizable CG water 85 may address some concerns surrounding CG force fields, the speed-up from atomistic is around three times less than classical MARTINI, thus rendering its use prohibitive in the case of the simulations performed here. It should also be noted that the use of polarizable water, did not improve the protein-lipid interaction in the case of aquaporin and GM3 in a recent study 79 .
In terms of convergence we performed tests based on protein MSDs (see SI Figure S4), which revealed that the protein MSD per frame reaches a plateau once large clusters (involving at least half of all channels in a membrane) formed, suggesting that our simulation lengths were sufficient to probe clustering and the resultant diffusional behavior of both channels and lipids.

Conclusions
Here we have presented insights into the nanoscale dynamic organization of cell membranes obtained from large-scale coarse-grained simulations of multiple Kir channels in model membranes at physiological degrees of protein crowding. Clustering of the channel proteins was modulated by the compositional complexity of the lipid bilayer, which altered not only the dynamics but also the spatiotemporal patterns of cluster formation.
Diffusion coefficients for both protein and lipid were lowest in the most complex lipid mixtures. Removing cholesterol caused the largest increase in diffusion coefficients, as might be anticipated 57,58 . Anomalous diffusion was observed for all proteins and lipids, and the anomalous diffusion of those lipids, which interacted strongly with the protein, was most affected by channel crowding; protein clustering slowed and perturbed diffusion of the interacting lipids (PIP 2 and GM3). Protein crowding and clustering also altered the mechanical properties of the bilayer (for instance, bilayer fluctuations) in a PIP 2 dependent fashion.
These simulations provide us with a glimpse of the complexities inherent in the interplay between the dynamic organization of Kir channels and of their allosteric activator lipid PIP 2 86 in cell membranes. As we have demonstrated, protein clustering can slow and cluster lipid modulators and form dynamic compartments in the membrane, which in turn may play an important role in the kinetics of channel activation 68 . More generally, our results imply that in order to understand the biological function of a cell membrane one must fully characterize the dynamic organization of the proteins and lipids on a nano to microscale.

Methods
System set-up and simulations performed. Simulations performed are shown in Table 1 and Table S1.
The Kir2.2 channel atomistic structure (PDB ID: 3SPI 32 ); was converted to a coarse-grained representation using the martinize.py script (version 2.4) downloaded from the MARTINI website (http://md.chem.rug.nl/index.php/ tools2/proteins-and-bilayers) with the MARTINI2.2 force field 87,88 and the ElNeDyn elastic network model 89 . A single CG model of Kir2.2 was embedded into a phosphatidyl choline (PC) bilayer by allowing lipids to self-assemble around the protein in a short simulation 90 . The PC membrane was subsequently converted into a membrane of the desired lipid composition using an in-house exchange lipid methodology 54 . The Kir2.2 channel in a model plasma membrane (PM) was then equilibrated for 100 ns prior to being tessellated, using the GROMACS tool genconf, to form 12 × 12, 6 × 6 and 3 × 3 grids. The resulting systems contained either 144, 36 or 9 Kir2.2 channels, ~55,000 (12 × 12 and 6 × 6 systems) or ~3,500 (3 × 3 systems) lipids, and were solvated using the standard MARTINI water model and neutralized to a 0.15 M NaCl concentration (see Table 1 for details).
Lipid bilayers were modeled with a mixture of the most abundant lipids present in mammalian cell plasma membrane 2-4 , with lipids distributed asymmetrically between the inner and outer leaflets, as previously described 54,55 . In the PM model, which is the most complex presented here, the outer leaflet contained phosphatidyl choline (PC)/phosphatidyl ethanolamine (PE)/sphingolipid (Sph)/ganglioside (GM3)/cholesterol (Chol) lipids in a ratio of 40:10:15:10:25, while the inner leaflet contained PC/PE/phosphatidyl serine (PS)/phosphatidylinositol 4,5-bisphosphate (PIP 2 )/Chol in a ratio of 10:40:15:10:25. PC, PE, and PS were modeled with 1-palmitoyl-2-oleoyl (i.e., PO) lipid tails. PIP 2 was removed from the PM model by substitution with PE, as these two species clustered in areas of higher membrane curvature in previous PM simulations 54 . In simulations where lipid complexity was further reduced, the proportions of the PM simulation were maintained for all remaining lipid species (Table 1 and Table S1). Sph and GM3 were modeled with a monounsaturated ceramide tail; PIP 2 was modeled using fully saturated tails.
Simulation parameters. All simulations were performed using GROMACS 4.6 (www.gromacs.org) and the standard MARTINI protocol. Periodic boundary conditions were applied, and a time step of 20 fs was used in all simulations. The temperature was maintained at 323 K using a Berendsen thermostat 91 , and the pressure at 1 bar using a Berendsen barostat, except in the case of the PCnoprot simulation, where a velocity rescale thermostat 92 and Parinello-Rahman barostat 93 were used. For both the temperature and pressure, a coupling constant of 4 ps was used for all simulations containing protein and 1 ps for the PMnoprot and PCnoprot simulations. In all simulations, the reaction field coulomb type was used with a switching function from 0.0 to 1.2 nm, and the van der Waals interactions were cutoff at 1.2 nm with a switching function applied from 0.9 nm. The LINCS algorithm 94 was used to constrain covalent bonds to their equilibrium values. All systems containing protein were simulated for between 20 and 50 μs, whereas systems without proteins were performed for 10 μs (PMnoprot) or 5 μs (PCnoprot); see Table 1.
Simulation analysis. Radial distribution functions of lipid headgroups around the channels were calculated using the GROMACS tool, g_rdf. Clustering analysis was performed using in-house Python scripts, making use of the NumPy 95 , MDAnalysis 96,97 and NetworkX 98 Python libraries. In the clustering algorithm, proteins were considered to be interacting when the centroids of the cytoplasmic domains (defined as residues 1-45 and 141-377) were within 8.2 nm of one another. Protein-protein interactions were identified using in-house clustering scripts, again making use of the NumPy and MDAnalysis modules. Residues of neighboring proteins were considered to be interacting when residue centroids were within 0.7 nm of one another. Phosphate densities were calculated using the g_density tool from GROMACS.
Diffusion analysis was carried out on trajectories with center of mass motion of the protein and lipids removed, and with the first 50 ns removed. Although the PM simulation was run for 50 µs, only the first 20 µs was used for diffusion analysis, in order to be comparable with analysis of the NoPIP2, PC, PMsparse and small simulations. Time-and ensemble-average mean squared displacement (MSD) for each species was extracted using a time step of 1 ns. Diffusion coefficients (D α ) were obtained by fitting the equation MSD = 4D α Δt α to graphs of MSD vs Δt, with Δt in the range 1 ns to 2 µs. Error bars show the standard deviation for the nonlinear fitting procedure. Scripts to perform this analysis were from https://github.com/tylerjereddy/diffusion_analysis_MD_simulations 99 . Alpha values were obtained by fitting a straight line at fixed intervals to the log-log plots of MSD vs Δt.
Power spectra of membrane undulations were estimated as detailed in 23 . The surface of the bilayer was calculated using lipid phosphate beads projected onto a grid of size 0.5 nm. The Fourier transform of the bilayer surface was calculated using FFTW routines in the SciPy Python library 100 , and 1D power spectra obtained by performing radial averaging. The code can be obtained from GitHub (https://github.com/philipwfowler/ calculate-bilayer-power-spectrum).