Lipid-Protein Interactions Are Unique Fingerprints for Membrane Proteins.

Cell membranes contain hundreds of different proteins and lipids in an asymmetric arrangement. Our current understanding of the detailed organization of cell membranes remains rather elusive, because of the challenge to study fluctuating nanoscale assemblies of lipids and proteins with the required spatiotemporal resolution. Here, we use molecular dynamics simulations to characterize the lipid environment of 10 different membrane proteins. To provide a realistic lipid environment, the proteins are embedded in a model plasma membrane, where more than 60 lipid species are represented, asymmetrically distributed between the leaflets. The simulations detail how each protein modulates its local lipid environment in a unique way, through enrichment or depletion of specific lipid components, resulting in thickness and curvature gradients. Our results provide a molecular glimpse of the complexity of lipid-protein interactions, with potentially far-reaching implications for our understanding of the overall organization of real cell membranes.

Introduction

Typical cell membranes are composed of hundreds of different lipid types that are asymmetrically distributed between the two leaflets. Embedded are many different membrane proteins, covering an estimated membrane area as large as 30% at a lipid/protein ratio of about 50–100.[1] The variety in cell membrane components gives rise to complex lipid–protein interplay.[2,3] Lipids do not simply provide the matrix where proteins are embedded but can actively participate in the regulation of protein activity, trafficking, and localization.[3] Proteins can either bind lipids specifically, where a clear binding site for a given lipid can be identified, or nonspecifically, where lipids act as a medium, and physical properties like thickness, fluidity, or curvature regulate protein function.[4,5]

The characterization of lipid–protein interactions is a key factor in our understanding of the organizational principles of cell membranes. Several experimental techniques are available to probe these interactions.[6] X-ray crystallography and electron crystallization can be used to identify lipids strongly bound to proteins as these lipids have to survive the crystallization process.[7,8] Lipid binding to membrane proteins can also be studied using fluorescence methods[9] or by mass-spectrometry on isolated lipid–protein complexes.[10] The recent development of extraction of membrane proteins from their native environment using nanodiscs is very promising in this regard.[11,12] Nevertheless, these techniques mainly capture strong interactions, and although some are quantitative, they do not give high spatial resolution.

Computational approaches on the other hand, such as molecular dynamics (MD) simulations, can provide such details and have been extensively used to study lipid–protein interplay.[13−16] In particular the use of coarse-grain (CG) models allows simulation of reversible binding and unbinding events and identify both strong and weakly binding lipids.[17,18] So far, most computational studies have been restricted to model membrane systems with a few lipid types. The recent modeling of a complex plasma membrane mixture containing more than 60 different lipids, however, has opened the way to probe lipid–protein interplay in a more realistic membrane environment.[19−22] Here, we extend this work by analyzing the lipids around 10 different classes of plasma membrane proteins. We find that each protein forms its own unique lipid shell, which gives rise to a complex and nonuniform perturbation of local membrane properties (“fingerprint”). The results show a rich variety of lipid–protein interactions and protein effects on membrane properties, emphasizing the importance of not just tightly bound lipids but the overall structure of the lipid–protein matrix.

Results and Discussion

Our simulation setup consists of a membrane patch containing around 6000 lipids, represented by the CG Martini force field.[23] The system measures ca. 42 × 42 nm in the lateral dimensions, and contains 63 different lipid types distributed asymmetrically between the leaflets. The membrane composition is based on previous work,[20] and models a prototypical plasma membrane. The outer leaflet contains most of the PC and SM lipids and all of the gangliosides (GM1 and GM3), has a higher level of saturation of the tails, and is slightly enriched in cholesterol. The inner leaflet contains all anionic lipids (e.g., PG, PIPs) and most of the PE, and has a higher level of polyunsaturation. Full details of the composition are given in the Supporting Information.

As shown in Figure A for aquaporin 1 (AQP1), we embedded four copies of the same membrane protein inside our plasma membrane model. The four proteins are weakly restrained at their initial position, at a distance of ca. 20 nm from each other, providing a computationally efficient way to increase statistics on lipid distribution around the proteins and to obtain an additional estimate of statistical errors independent of time correlations. Based on 30 μs long simulations, we characterize the distinctive nature of the lipid environment surrounding each protein, using a distance cutoff criterion of 0.7 nm (Figure B). Averaged over the last 5 μs of simulation time, we then compare the composition of this lipid shell around the four copies of the protein to the bulk plasma membrane composition, expressed as the relative depletion–enrichment (D–E) index for different categories of lipids (Figure C). The same analysis was repeated for lipid shells of 1.4 and 2.l nm around the AQP1 molecules (Table S1).

Figure 1

Unique lipid environments for different membrane proteins. (A) Simulation setup, consisting of a plasma membrane model containing 63 different lipid types with four membrane proteins embedded. (B) View of the local lipid environment around AQP1, displaying lipids within a distance cutoff of 0.7 nm from the protein surface. (C) Lipid depletion–enrichment (D–E) index in the case of AQP1, obtained from the last 5 μs of a 30 μs long simulation, and averaged over the four AQP1 molecules (error bars indicate standard deviation). The D–E index is computed by dividing the lipid composition of the first 0.7 nm shell by the bulk membrane composition. Values larger than 1 indicate enrichment of a given lipid group, while values smaller than 1 indicate depletion. (D) D–E index matrix, with average depletion (blue)/enrichment (red) for 10 different membrane proteins (the calculation for additional distance cutoffs and corresponding standard deviations are shown in Table S1). The COX-1 D–E index values for the negatively charged lipids of the lower leaflet have been omitted because they are difficult to interpret given the partial insertion of the protein only in the upper leaflet (see the note in Table S1). Lipid classes considered are phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), phosphatidic acid (PA), diacyl-glycerol (DG), lyso-PC (LPC), sphingomyelin (SM), ceramide (CER), phosphatidylinositol (PI), phosphatidylinositol-(bi, tri)phosphate (PIPs), ganglioside (GM), cholesterol (CHOL), polyunsaturated (PU), fully saturated (FS), and others.

To obtain a general view of the diversity of lipid shells around membrane proteins, we performed this simulation protocol for 10 diverse plasma membrane protein families (Figure S1). The proteins considered, besides AQP1, are prostaglandin H2 synthase (COX1), dopamine transporter (DAT), epidermal growth factor (EGFR), AMPA-sensitive glutamate receptor (GluA2), glucose transporter (GLUT1), voltage-dependent Shaker potassium channel 1.2 (Kv1.2), sodium, potassium pump (Na,K-ATPase), δ-opioid receptor (δ-OPR), and P-glycoprotein (P-gp). These proteins include transporters, channels, enzymes, and receptors, and represent different quaternary structures and sizes, as well as monotopic membrane proteins (Figure S1). As a result of the differences in sizes between the proteins, the amount of lipids found in the nearest 0.7 nm shell varies significantly, from ca. 32 as in EGFR (whose transmembrane domain consists of only two helices) to ca. 78 for AQP1 and ca. 95 for Kv1.2 (which are tetrameric proteins). To compare the lipid shells of the different proteins irrespective of protein size, we compute the lipid D–E index for each of the 10 systems. The results are presented in Figure D and Table S1 for all 10 proteins and the three distance cutoffs used to define the lipids shells.

The D–E index values, evaluated in terms of statistical significance, reveal some interesting generic trends (Supporting Information, and Tables S1–S9). Within the first lipid shell, approximated using the distance criterion cutoff of 0.7 nm, we observe a depletion of PC and SM lipids across the whole spectrum of membrane proteins (Tables S1 and S6). Additional control simulations performed on the AQP1 system to test for convergence of the results in terms of simulation conditions and length confirmed this trend (Tables S2–S5 and S7–S9). A class of lipids for which we detect a significant enrichment is represented by the gangliosides (GM lipids), present only in the upper leaflet. Gangliosides are key signaling lipids and are known to bind to proteins and regulate their sorting behavior and activity.[24−30] In the simulations, GM enrichment extends at least up to 2 nm from the proteins (Tables S1–S9). In the lower leaflet, enrichment of PIP lipids is common to all the simulation systems (with the exception of COX1, a monotopic protein adsorbed at the outer plasma membrane leaflet). PIPs, too, are important signaling lipids and are found bound to a large variety of proteins,[31−34] and implicated in controlling protein–protein interactions.[35] Although the enrichment of GM and PIP lipids depends somewhat on the water model used in the simulations (Tables S2), it is detected in all the systems and for different simulation setups (Tables S1–S9).

Despite these common trends, the most striking conclusion from our D–E index analysis is the large differences found between the individual protein classes. Each protein has a unique composition of the first lipid shell. For instance, the enrichment of PIP lipids can be accompanied by a certain degree of enrichment of another class of negatively charged lipids, the PI lipids, with the exception of GluA2 and δ-OPR, while the results for PS and PA lipids, also negatively charged, show larger differences across proteins. For DAG molecules we observe a significant enrichment near the proteins, with the exception of Kv1.2 and P-gp.

When lipids are grouped based on degree of unsaturation, in fully saturated (FS) and polyunsaturated (PU; see the Supporting Information for full details on the definition) classes, we observe a significant enrichment of PU lipids for several proteins, including DAT, GluA2, GLUT1, Kv1.2, Na,K-ATPase, and P-gp, but not for all proteins. The results for FS lipids show a less clear trend, with a D–E index value often very close to 1, and no significant enrichment or depletion within the first lipid shell. We also observe that the enrichment or depletion of a given class can be correlated with membrane composition: we have simulated the AQP1 system in the absence of GM lipids, which resulted in a more significant enrichment of PU lipids near the proteins compared to the simulations with GM lipids (Tables S2–S3 and S7–S8) but similar results for other classes. However, a similar pattern was not detected for the Na,K-ATPase, simulated for an additional 20 μs with no GM lipids after the 30 μs long simulation in the presence of GM (Tables S1, S5, S6, and S9). For this system, in fact, the enrichment of PU lipids persisted with or without GM lipids, highlighting the different behavior of each protein in the plasma membrane mixture.

The variation in cholesterol concentration in the lipid shells is also remarkable. Note that cholesterol is present at an overall probability of 30%, implying a local concentration of 35–40% in the case of Kv1.2, AQP1, and δ-OPR (D–E index 1.2–1.3), as opposed to a concentration of around 25% in the case of DAT (D–E index 0.8). For DAT, depletion of cholesterol is accompanied by an enrichment in PU lipids and GM lipids, many of which are fully saturated, while for Kv1.2 we detect enrichment in cholesterol, PU, and GM within the same shell. Despite the common notion that cholesterol prefers to interact with fully saturated lipids and in fact induces phase separation,[36−38] the coexistence of sites to specifically accommodate cholesterol and PU lipid molecules on the protein surface has also been previously shown for rhodopsin.[39,40] Together, our data show how proteins can have a strong impact on the lateral organization of lipids in their immediate surrounding.

However, the impact of protein on membrane organization is even more complex. The preceding analysis only revealed information on the average composition of the lipid shells up to 2.1 nm from the proteins. Considering that membrane proteins are not simple homogeneous cylindrical inclusions, one can anticipate that the distribution of lipid species around the protein is nonuniform.

To investigate this, we analyzed the local density fluctuations around each of the proteins by computing 2D lateral density maps. These maps were obtained by counting the presence of lipid beads in grid cells with 0.3 nm spacing and averaging over the last 5 μs of the simulations. Given the complex composition of our membrane, we restricted our analysis here to three major lipid classes, namely, PU lipids, FS lipids, and cholesterol (CHOL). For each class, depletion or enrichment are here shown with respect to their corresponding average density calculated in a given leaflet. The results for AQP1 are shown as an example (Figure A), but similar density fluctuation patterns are seen around all proteins (Figure S2).

Figure 2

Membrane protein fingerprints. (A) Two-dimensional lateral density maps, showing local density fluctuations around AQP1 in upper (top row) and lower (bottom row) leaflets, grouped according to lipid classes: polyunsaturated (PU) lipids, fully saturated (FS) lipids, and cholesterol. Major observations are indicated by arrows, see text for details: I, nonspecific binding; II, nonuniform distribution; III, leaflet asymmetry; IV, specific binding; V, membrane fluctuations. (B) Nonuniform variations in local membrane properties around AQP1: thickness, mean curvature, and Gaussian curvature for upper (top row) and lower (bottom row) leaflets.

We observe a rich spectrum of possible features of protein–lipid interactions, indicated by arrows I–V in Figure A. One important feature is the broad distribution of FS and PU lipids, in both leaflets, near many proteins, including AQP1, DAT, EGFR, GluA2, GLUT1, Na,K-ATPase, and δ-OPR (Figure A, arrow I, and Figure S2). Other clear features are that the observed density fluctuations are strongly inhomogeneous around the protein (arrow II, Figure A), as well as leaflet-specific as they depend on protein structure and lipid composition, both of which are asymmetric (arrow III, Figure A).

In the upper leaflet, for example, we notice regions of strong FS enrichment in contact with the proteins. Such regions are usually more localized than PU enriched regions and are often coupled with smaller, yet still highly localized, FS enriched regions in the lower leaflet. This behavior can be seen for AQP1, DAT, EGFR, δ-OPR, and even for the monotopic COX1, which is only partially embedded in the upper leaflet (Figure A and Figure S2). The size and shape of FS lipid regions differ from protein to protein, and in many cases create a discontinuous ring around the transmembrane domains, as in the case of AQP1, GluA2, and Kv1.2, which are homotetramers (Figure A and Figure S2).

In the lower leaflet, PU enrichment is often observed near the proteins, and particularly noticeable around, for example, DAT, GluA2, GLUT1, and Na,K-ATPase (Figure A and Figure S2).

Some proteins clearly induce a sharp partitioning of the different lipid classes. This is the case for GluA2 and PU lipids in the lower leaflet, and P-gp, where in the upper leaflet we observe a clear distinction between the side of the transmembrane domains in contact with PU lipids and the side in contact with FS lipids, while in the lower leaflet there is an obvious preference for PU lipids (Figure S2). Kv1.2 is another striking example of how the same lipid class (PU) can be asymmetrically distributed between leaflets, and symmetrically distributed around the protein within the same leaflet, due to the homotetrameric nature of the channel and possibly linked to a more specific type of interaction (Figure S2).

Monotopic proteins are also capable of inducing a clear separation in the distribution of lipid classes, even in the leaflet they are not directly bound to. For COX1, for example, in the lower leaflet the enrichment of PU lipids is aligned with the depletion of FS lipids and CHOL underneath the protein, partially embedded only in the upper leaflet.

In addition, to such a nonspecific, broad distribution of lipid classes, we detect an example of more specific interactions, especially for the CHOL class (arrow IV), where the term specific in this case is used to indicate a localized region of enrichment in close proximity with the protein, possibly as a consequence of longer contact times between the protein and the cholesterol molecules due to the presence of preferred cholesterol interaction sites. Cholesterol is the most abundant single component of the plasma membrane mixture, and thus associated with a more even distribution compared to the PU and FS classes. However, the corresponding 2D density maps suggest sites of possibly specific cholesterol interactions. AQP1, for example, clearly shows localized enrichment of cholesterol at the interface between monomers, and similar features are also observed in Kv1.2, DAT, GLUT1, Na,K-ATPase, and δ-OPR. The D–E index analysis shows, for some proteins, enrichment of cholesterol and PU lipids within the same shell (Figure and Table S1). The 2D density maps highlight how such regions, when present, do not substantially overlap in space (Figure and Figure S2).

For proteins such as GLUT1, Na,K-ATPase, and DAT the presence of specific motifs or sites for cholesterol binding is reported in the literature.[41−43] Here, we focus on DAT, as an example of protein for which our simulations retrieve the experimental cholesterol binding site, and on Kv1.2, for which we predict specific cholesterol interactions. In the crystal structure of DAT a cholesterol molecule is bound in a groove lined by residues of transmembrane helix 1, 5, and 7, closer to the cytoplasmic side of the membrane.[44] In our simulation, the presence of cholesterol is clearly visible for all four DAT molecules, in the same groove that coordinates the cholesterol molecule in the crystal structure (Figure S3A). The simulation, moreover, reveals the dynamic component of cholesterol interactions at this site, as multiple binding modes are detected (Figure S3A), in line with what has been already described for the corresponding cholesterol site in the serotonin transporter.[45]

For Kv1.2, we observe cholesterol molecules in a site located at the interface between two monomers and protruding toward the central pore of the channel (residues L328, L331, I402, P405, and V406 of one monomer, and I396, A397, and L400 of the neighbor one) (Figure S3B). For a different type of potassium channel, the inward-rectified Kir2.1 potassium channel, specific residues able to modulate the sensitivity of the channel to cholesterol are located at the interface between monomers.[46] The sequence alignment of several Kir channels with Kv1.2 shows that the residues’ cholesterol-binding residues here identified for Kv1.2 correspond to some of the residues that regulate cholesterol sensitivity of Kir2 channels (in particular, Kir2.1-I166 and Kir2.1-I175 correspond to Kv1.2-I396 and Kv1.2-P405, respectively).[43] Thus, our CG simulations further strengthen the hypothesis that such cholesterol interaction sites might have a functional role that spans across different types of potassium channels.[46] These results also suggest that our computational approach might help in identifying protein sites primed to interact with specific lipids, even when experimental data are not available.

Finally, significant lipid density fluctuations can be seen extending to the rest of the membrane, away from the protein surface (arrow V, Figure A). Similar density fluctuations were also observed in the pure plasma membrane models,[20,47] but it is likely that the proteins affect the extent and nature of these. However, given the currently used setup with four copies of each protein being present in the simulation cell, analysis of the long-range effects of these density fluctuations is problematic and left for future studies.

One possible reason for proteins to accumulate certain classes of lipids around them is to ensure a proper embedding of their hydrophobic transmembrane domains.[48] Thus, one would expect variations in local thickness profiles and curvature gradients, as has been predicted by theory,[49] and observed in many previous MD simulations of proteins in simpler membranes.[21,50−54] The 2D membrane thickness profile around AQP1 is shown in Figure B, alongside the mean and Gaussian curvature fields. Indeed, we observe strong perturbations in all of these properties, with similar perturbations around all 10 proteins (Figures S4–S6). Again, we retrieve the same type of inhomogeneity as for the density maps, with perturbations strongly varying along the surface of the protein and depending on the leaflet. Note, for instance, the strong alteration in both mean and Gaussian curvature around the rim of AQP1 (Figure B).

MD simulations have been used extensively to study how lipids might regulate protein function, and recent advances in high-performance computing and in the development of CG models have allowed for a higher degree of complexity in simulated systems.[13,20−23,55−61] In the Martini model used in this study,[62−64] a compromise between computational efficiency and chemical details of the system is achieved by grouping an average of four heavy atoms into a larger particle or bead: several bead types and subtypes exist, to account for the chemical diversity of the systems to model. In the case of proteins, the lower resolution provided by the Martini model requires the use of an elastic network to preserve their tertiary structure,[65] and as a consequence, the force field does not allow the study of conformational transitions of proteins, which could be due, for example, to lipid binding. Moreover, the tendency of the force field to overstabilize protein–protein interactions and form protein–protein aggregates that do not dissociate over long time-scales has recently been reported.[21,66] Here, with these caveats in mind, we used the Martini model to describe lipid–protein interactions in a plasma membrane mixture over 30 μs long simulations, which are not yet feasible at an atomistic level of detail due to the complexity of the systems. To account for reproducibility of the lipid distribution around a given protein type, our approach considered four molecules of the same protein, restrained to avoid the formation of protein–protein aggregates that cannot be sampled efficiently. We have tested different simulation setups and lengths, analyzed the different profiles in lipid enrichment near the proteins (Tables S1–S9), and obtained qualitatively similar results.

Our simulations reveal a hitherto unappreciated level of complexity of protein–lipid interplay in the asymmetric environment of a realistic plasma membrane. It should be considered, however, that the details of the lipid compositions of the first lipid shells near the proteins found in our study are, to some extent, sensitive to the details of the force field used, although most of the features we observed are plausible based on the results of other studies. For instance, GM lipids are known to contribute to the lateral organization of membranes, due to their large head groups, their saturated tails, and high phase transition temperature.[67,68] We observed small differences in the extent of GM aggregation around the proteins depending on the water model used as also pointed out in a previous work.[24] However, the effects of the presence of GM lipids in membrane models have been studied using both atomistic and CG simulations (and reviewed, for example, by Manna et al.),[69] with atomistic simulations highlighting the tendency of GM molecules to interact with each other, and CG simulations showing their ability to segregate in nanodomains in model membranes.[25,60,70] Experimental work on protein–GM lipid interactions includes, among others, raft-associated glycophosphatidylinositol-anchored proteins (where protein enrichment is colocalized with GM enrichment),[71] and growth factors.[72] For EGFR in particular,[30] CG simulations were used to characterized in more detail the interactions with GM lipids,[27] and evaluate the strength of such interactions.[73] Additional details on GM–protein interactions are discussed in the Supporting Information.

The clustering of negatively charged lipids around membrane proteins is affected by the electrostatic interactions’ treatment in the simulations. However, their enrichment, despite differences in magnitude for PIP, for instance, is present for all the proteins, and has also been reported by others. Several works, for example, focus on the interactions and redistribution of cardiolipins around membrane proteins, providing results in agreement with available crystal structures (see, for instance, refs (74 and 75)), while others address the interactions between negatively charged lipids and ion channels,[17] tyrosine kinase receptors,[17,76,77] or GPCRs.[21] Multiscale simulation approaches were used to characterize the interactions between the juxtamembrane domain of EGFR and PIP lipids, in bilayers with different ratios of POPC, POPS, and PIP2 lipids,[76] and also reported PIP lipid redistribution and interactions with the ion channel Kir2.2 in agreement with experimental data.[78] Additional details on PIP lipid–protein interactions are discussed in the Supporting Information.

The lateral organization of PU lipids has also been addressed by simulations. For the membrane model used in this work and a similar membrane model, PU lipids form clusters of the size of few lipids only and without inducing any phase separation.[20,47] This excludes that the regions enriched in PU lipids observed here near the membrane proteins are biased by the force field parameters. And indeed, for simpler mixtures the Martini model reproduces experimental findings showing that PU lipids promote the coexistence of liquid ordered and disordered phases by accumulating and fluidizing the liquid disordered phase[79,80] and shows free energies of transfer for polyunsaturated DLiPC lipids between different environments that are comparable between Martini and atomistic simulations,[81] although with different entropy and enthalpy contributions as expected from a coarse-grained model.[23] The presence of regions enriched in PU lipids or polyunsaturated fatty acids (PUFAs) near membrane proteins has been observed in atomistic and CG simulation studies before.[39,40,82−85] Of note are atomistic simulations of open and closed states of the Shaker Kv channel, which have shown accumulation of PUFAs on the upper channel surface, similarly to what we observe in our CG simulations of the Kv1.2 channel.[85]

Despite the simplified representation of the system, the Martini model has been successfully applied to elucidate lipid binding sites around several examples of membrane proteins, including recovering of known experimental binding sites as well as predicting novel ones.[17,22,74,76,78,86−90] Here, we find that the distributions of lipid types around each protein, and near different parts of each protein, are distinct, providing a unique environment or “lipid fingerprint” for each protein. These unique lipid shells, in turn, affect local membrane properties such as thickness and curvature gradients in a highly nonuniform way. We expect also other membrane properties, such as the lateral pressure profiles or local compressibility as well as dynamics, to be strongly anisotropic. The lipid fingerprints offer a large range of opportunities for cells to steer and control membrane organization. For instance, specific binding sites could act as glue, stabilizing protein–protein contacts, or cause blockage of interactions with other proteins. Similarly, merging of nonspecific binding sites could recruit different proteins into close proximity forming nanosized domains. The protein-induced curvature gradients may provide long-range effects, driving either protein oligomerization or repulsion. The above sorting mechanisms are already known, and used in a variety of mean field models to predict large-scale membrane organization,[91−94] but our data suggest such models are too simplistic. The interactions between membrane proteins will strongly depend on the protein pair considered, and be highly nonuniform, i.e., depending on their relative orientation. For instance, proteins may attract each other in certain orientations, as shown by studies on the ATP-synthase,[95,96] but repel each other in others. Furthermore, through control of the lipid fingerprint, these interactions will depend very sensitively on the overall composition of the membrane. Moreover, in the crowded environment of typical plasma membranes, another level of complexity is added through the competition for binding sites between the proteins. Taken together, we observe a new level of complexity in protein–lipid interplay, with profound implications for the lateral and dynamic organization of cell membranes. Furthermore, better understanding of the nature of the local environment might help with designing and understanding membrane mimetics for functional studies of membrane proteins.[97] We therefore expect that detailed simulations of increasingly realistic cell envelopes, together with improved characterization of lipid shells with experimental methods, are needed to fully appreciate the role of lipid fingerprints in cellular processes.

Methods

The detailed methods are described in the Supporting Information. In summary, we embedded 10 membrane proteins in a previously characterized model of the plasma membrane.[20] The starting structures of the 10 membrane proteins simulated in this study were taken from the Protein Data Bank or obtained from the corresponding publication: aquaporin-1 (AQP1, PDB ID 1J4N);[98] prostaglandin H2 synthase (COX1, PDB ID 1Q4G);[99] the dopamine transporter (DAT, PDB ID 4M48);[44] the epidermal growth factor receptor (EGFR);[77] AMPA-sensitive glutamate receptor 2 (GluA2, PDB ID 3KG2);[100] glucose transporter 1 (GluT1, PDB ID 4PYP);[101] voltage-dependent Shaker potassium channel 1.2 (Kv1.2, PDB ID 3LUT,[102] residues 32 to 4421 for each monomer); sodium, potassium pump (Na,K-ATPase, PDB ID 4HYT);[103] δ-opioid receptor (δ-OPR, PDB ID 4N6H);[104] and P-glycoprotein (P-gp, PDB ID 4M1M).[105] In each system, four copies of each protein were included and positioned at a distance of ca. 20 nm from each other. Proteins were simulated using standard Martini protocols with minor variations between systems to accommodate system-specific issues (Supporting Information). The following lipid classes were included: cholesterol (CHOL), in both leaflets; charged lipids phosphatidylserine (PS), phosphatidic acid (PA), phosphatidylinositol (PI), and the PI-phosphate, PI-bisphosphate, and PI-trisphosphate (PIPs) placed in the inner leaflet; and ganglioside (GM) in the outer leaflet. The zwitterionic phosphatidylcholine (PC), phosphatidylethanolamine (PE), and sphingomyelin (SM) lipids were placed in both leaflets, with PC and SM primarily in the outer leaflet and PE in the inner leaflet. Ceramide (CER), diacylglycerol (DAG), and lysophosphatidylcholine (LPC) lipids were also included, with all the LPC in the inner leaflet, and CER and DAG primarily in the outer leaflet. The details of the Martini lipids used in this study can be found on the Martini Lipidome webpage (http://www.cgmartini.nl/index.php/force-field-parameters/lipids) and are described by Ingolfsson et al., and Wassenaar et al.[20,106] The exact lipid composition of each system is given in the Supporting Information. The systems are ca. 42 × 42 nm in the membrane plane (x and y), including 4 proteins and ca. 6000 lipids.

Production runs were performed in the presence of weak position restraints applied to the protein backbone beads, with a force constant of 1 kJ mol–1 nm–2, preventing proteins from associating with each other. Each of the systems has been simulated for 30 μs, which turned out to be adequate to obtain convergence of major lipid components in the lipid shells around the individual copies of the proteins (Supporting Information). Additional control simulations were performed in the AQP1 system, in order to test the effects of simulation length, position restraints on the proteins, lipid composition, and water model on the results of lipid composition near the proteins (Supporting Information).

Simulations were performed using the GROMACS simulation package version 4.6.3,[107] with the Martini v2.2 force field parameters,[62,63] and standard simulation settings.[108] Additional details are provided in Supporting Information. All the analyses were performed on the last 5 μs of each simulation system.