Capturing Biologically Complex Tissue-Specific Membranes at Different Levels of Compositional Complexity.

Plasma membranes (PMs) contain hundreds of different lipid species that contribute differently to overall bilayer properties. By modulation of these properties, membrane protein function can be affected. Furthermore, inhomogeneous lipid mixing and domains of lipid enrichment/depletion can sort proteins and provide optimal local environments. Recent coarse-grained (CG) Martini molecular dynamics efforts have provided glimpses into lipid organization of different PMs: an "Average" and a "Brain" PM. Their high complexity and large size require long simulations (∼80 μs) for proper sampling. Thus, these simulations are computationally taxing. This level of complexity is beyond the possibilities of all-atom simulations, raising the question-what complexity is needed for "realistic" bilayer properties? We constructed CG Martini PM models of varying complexity (63 down to 8 different lipids). Lipid tail saturations and headgroup combinations were kept as consistent as possible for the "tissues'" (Average/Brain) at three levels of compositional complexity. For each system, we analyzed membrane properties to evaluate which features can be retained at lower complexity and validate eight-component bilayers that can act as reliable mimetics for Average or Brain PMs. Systems of reduced complexity deliver a more robust and malleable tool for computational membrane studies and allow for equivalent all-atom simulations and experiments.

Introduction

Membrane proteins constitute a disproportionately large fraction of all drug targets,[1] as they are often the gate-keepers of cellular signaling and regulation. As such, detailed understanding of membrane protein function is critical for the design and development of future drug candidates. Working experimentally with membrane proteins is extremely difficult. In order to retain the correct structure and function of the protein during ex vivo experiments, one must reconstitute the protein either into an artificial membrane environment or into a suitable mimetic (such as a nanolipoprotein particle[2−5] or micelle[6]).

For the most part, the function of these in vitro platforms for membrane protein study is to provide a hydrophobic environment that mimics the center of a lipid bilayer to maintain the secondary structure elements that are specific to integral membrane proteins (such as long α-helices that either are amphipathic or have many hydrophobic residues on their surface). Even many early membrane protein simulation studies[7,8] did not use phospholipids but rather utilized organic solvents such as octanol, as they were computationally less complex, were better parametrized, and allowed for faster relaxation and equilibration of the system—an essential necessity when computational resources were limited. Recent advancements in both computational and experimental methods for studying membrane proteins have revealed that the membrane is contributing to function and doing far more than simply acting as a “hydrophobic slab” in which to embed the protein and stop the secondary structure from collapsing.

The properties of the membrane and the properties of the embedded membrane protein are intrinsically linked. Membrane proteins are found to preferentially co-localize and interact with certain lipid types;[9−14] the local lipid environment rearranges based on properties such as lipid charge,[12,15−17] lipid tail order,[12,18,19] or hydrophobic thickness.[20−22] The notion of “lipid fingerprints” has recently been explored,[23] the idea that different proteins can have an identifiable pattern of lipids around them. Lipids can modulate the behavior of the protein,[9,10,12,15] and indeed, certain lipid types may even be a requirement for protein function.[9,15]

The complex nature of the membrane–membrane protein symbiosis becomes even more intricate when small chemical drug compounds are introduced. A new category of compounds has been defined as PAINS (Pan Assay Interference compoundS), or PAINS-like molecules, that were found to cause false-positives in many membrane protein assays by partitioning into the membrane and triggering a change in protein behavior through modulation of the membrane properties.[24−26] Thus, it is not just the protein but the nature of the membrane environment that has to be considered when studying membrane protein function. Consequently, the outcome of protein–drug interactions may be modulated by the local lipid environment and ultimately lead to different tissue-specific outcomes. For more reliable and beneficial computational contributions to study membrane protein function, there thus exists the need to use more realistic mixtures of lipids that can act as good mimetics for plasma membranes and that are even tissue- and/or disease-state-specific.

In attempts to more accurately represent biological plasma membranes, there has been a huge growth in the complexity of the lipid mixtures used for simulation. There are several recent reviews on computational modeling of complex membranes.[27−29] Work by the authors has also attempted to recreate the intricacy of different tissue types using ∼60 different lipid species.[30,31] However, in order to reach reasonable relaxation and equilibration of these massively complex lipid membranes, it required systems of ∼20,000 lipids in size and simulations 80 μs long. Both of these factors necessitated the use of the Martini[32,33] coarse-grained (CG) force field. To set up, run, and analyze these simulations necessitated considerable resources (both human and computational). There are many biological phenomena that require atomistic chemical detail in order to accurately simulate. Martini proteins need to have their secondary structure constrained, limiting the model’s use for studying membrane protein conformational changes. While Martini has been used extensively to study membrane domains,[34,35] it has been shown that the thermodynamic driving forces for cholesterol–phospholipid interactions are not captured with Martini.[36] Other limitations include a proper representation for water, limiting Martini’s applicability to study detailed interfacial properties that are crucial for phenomena such as pore formation[37] and the electrostatic potential. However, the length, size, and number of lipids in these large tissue-mimetic simulations mean that it is currently not feasible to implement an AA membrane simulation of 60 different components.

In this study, we demonstrate different levels of reduced complexity systems that maintain the general “tissue-specific” properties of each membrane. The reduced complexity of the system allows for smaller simulations that equilibrate faster. By validating against the “gold standard” highly complex 60-component membranes, we provide less complex mimetics that reproduce most of the membrane properties of their more complex counterparts. These specific environments are of importance for studying the behavior of different membrane proteins, as well as the different membranes themselves. Ultimately, the complexity of the membranes is reduced to such a level that it is feasible to convert the system from CG to AA representations.

Previously, we have shown that the complex (∼60 different lipid species) mixtures of lipids that represent different “tissue-like” membranes (an “Average” mammalian plasma membrane[30] and a representative “Brain” neuronal membrane[31]) have distinct properties that arise from their varying compositions. Teasing out the precise lipids responsible for the differing behaviors is challenging, especially as most properties are a cumulative effect from many different lipid types. However, there are some general membrane properties that can be measured and have observable differences between the two systems. As such, we decrease the complexity of our systems, while analyzing them to ensure that the general membrane properties are retained for each membrane type.

Our aim is to illustrate that lipid compositions of reasonable complexity can recapitulate the tissue-specific behaviors of plasma membranes, and thus provide robust membrane mimetics for future studies.

Methods

CG Simulation Setup

The bilayer builder insane(38) was used to construct the membranes. For each bilayer, the number of lipids in the inner and outer leaflets was adjusted based on an independent bilayer simulation using symmetrical compositions of either the inner or outer leaflet. This process was iterated with changes in the outer/inner leaflet distribution of cholesterol until the cholesterol distribution did not drift over time. This protocol is described in more detail in previous publications.[30,31,39] For a complete list of simulations and lipid compositions, see Table and the spreadsheet in the Supporting Information.

Table 1

Simulations Run

tissue type	lipid species	simulation name	no. of lipids	simulation time
Average (A)	63	A-63a,b	∼20,000	80 μs
	18	A-18	∼20,000	80 μs
	8	A-8	∼6000	40 μs
	8	A-8500	∼500	15 μs
	8	A-8AA	∼500	200 ns (atomistic)
Brain (B)	58	B-58b	∼20,000	80 μs
	16	B-16	∼20,000	80 μs
	8	B-8	∼6000	40 μs
	8	B-8500	∼500	15 μs
	8	B-8AA	∼500	200 ns (atomistic)

From ref (30).

From ref (31). All simulations are Martini CG unless specified otherwise.

CG Simulation Parameters

The CG simulations were performed using the Martini model[32,33] The lipid force fields used were described in Ingólfsson et al.[31] All of the lipid force fields can be found at the Martini portal www.cgmartini.nl. All of the simulations were run using either the GROMACS 4.6 or 5.1.4 simulation packages,[40,41] following the same setup described in Ingólfsson et al.[31] Briefly, a 20 fs time step was used for all production simulations with the standard Martini cutoffs, and the same “common” parameter set described in de Jong et al.[42] The simulations contained either ∼20,000 lipids or ∼6000 lipids with >15 CG waters per lipid (one CG water representing four water molecules), counterions, and 150 mM NaCl.

Undulations in the simulations were restricted using weak position restraints on selected lipids in the outer leaflet, as with previous studies.[30,31] The pressure and temperature were controlled using the Parrinello–Rahman barostat[43] (1 bar semi-isotropic pressure, with τp = 5.0 ps) and the velocity rescaling thermostat[44] (at 310 K, with τT = 1.0 ps). The larger membranes were simulated for 80 μs, while the smaller systems were simulated for 40 μs. Due to the smoother interaction potentials in Martini, the effective time of many processes is faster; i.e., for single component bilayers, lipid diffusion is roughly 4-fold faster at the Martini CG level.[32] This would correspond to 320 and 160 μs of effective time diffusion time in the large and smaller systems, but as this speedup is not universal, all reported times are unscaled.

Backmapping Protocol

Even the “smaller” CG systems of ∼6000 lipids (simulations A-8 and B-8, Table ) are unpractical for simulation using all-atom force fields. Thus, more feasible and appropriately smaller CG systems of both the A-8 and B-8 compositions were constructed and simulated using the same protocols as described above. These smaller systems were ∼500 lipids in size and were equilibrated for 15 μs (simulations A-8500 and B-8500, Table ). The final frames of these simulations were used as inputs to convert the CG system to AA resolution using the backmapping procedure described by Wassenaar et al.;[45] the initial CG to AA mapping is done by the backward.py script that places atoms based on defined bead mapping and mapping rules specified for each lipid type, and the initram.sh script uses the AA force field to run a number of minimization on MD steps to grow the initial structure to one that satisfies the AA force field. Due to the ∼4:1 mapping nature of Martini, multiple lengths of hydrocarbon tails are condensed into equivalent numbers of CG beads. The equivalent atomistic representations for saturated Martini tails of 3, 4, and 5 beads are C12:0–C14:0, C16:0–C18:0, and C20:0–C22:0, respectively (http://cgmartini.nl/index.php/force-field-parameters/lipids2/350-lipid-details). This means that a single CG tail can represent several potential tails in all-atom representations. Thus, for the CG lipid species in the A-8 and B-8 mixtures, the most appropriate AA lipid tails were chosen (Table S1). Missing Martini to CHARMM mapping files were constructed based on files from similar lipids and as specified in Wassenaar et al.[45] Note that rules were added to the mapping files to promote correct stereochemistry in initial structures but tests were not done to verify that those were maintained in all cases; therefore, these simulations should only be viewed as proof-of-principle.

All-Atom Simulation Parameters

The all-atom simulations used CHARMM36 lipids[46] and TIP3P water[47] and were simulated with GROMACS v5.1.4.[40] Simulations used a 2 fs time step. The Lennard-Jones interactions were cut off after 1.2 nm, with a switch-function from 1.0 to 1.2 nm. The particle mesh Ewald method[48,49] was used to calculate long-range electrostatic interactions. The pressure and temperature were controlled using the Parrinello–Rahman barostat[43] (1 bar semi-isotropic pressure, compressibility of 4.5 × 10–5, with τp = 12.0 ps) and the Nosé–Hoover thermostat[50] (at 310 K, with τT = 1.0 ps). Hydrogens were constrained with LINCS.[51]

Analysis Tools

For analysis of domains within the CG systems, the molecular configuration (3D coordinates) of the bilayer was converted to a surface representation using MemSurfer[52]—an open-source tool for characterizing and analyzing membrane surfaces. In particular, smooth surfaces were computed through optimization to best represent the lipid headgroups in the inner and outer leaflets of the membrane. The resulting surfaces were represented as periodic triangulated meshes, with lipid headgroups forming vertices of the meshes. Given these membrane surfaces, the probability density of cholesterol, γ(x), was estimated using the 2D kernel density estimation (KDE) approach with a Gaussian kernel, i.e.,where x represents the position of lipid headgroups and Nl and Nc are the number of all lipids and cholesterols in the leaflet, respectively. Here, the NC3 (or equivalent) CG beads were used to represent the positions of lipid headgroups, and the smoothing kernel σ = 3 nm was used for KDE. In the limit, γ(x) represents the probability of finding a cholesterol at any given point in the membrane and sums to 1. To obtain the cholesterol distribution function at any given position, the probability distribution was normalized, i.e., . The cholesterol distribution function represents the distribution of cholesterols in the domain and sums to the total number of cholesterols in the system. Hereafter, we take the cholesterol distribution function to represent the density of cholesterol in the system.

Although the density was estimated using 2D coordinates of lipid headgroups, a key reason to compute membrane surfaces is to explore the spatial variations of cholesterol distribution in order to identify domains. To this end, the topology of the cholesterol distribution fields (defined on triangular meshes) was analyzed. Connected regions above or below a certain value were taken to represent as domains, representing enrichment and depletion, respectively. Suitable thresholds to define such domains were identified through analysis. Given these thresholds of cholesterol distribution, characteristics of resulting domains were computed individually for each frame and then aggregated over all frames for each CG system. To efficiently explore the impact of different threshold choices as well as to compute statistics of interest, the topological analysis framework TALASS[53,54] was used to encode the statistics of all possible domains for all possible thresholds.

The reported average area per lipid (APL) for each lipid mixture was determined from the APL of flat symmetrical outer/inner systems for each composition. A force constant of 100 kJ mol–1 nm–2 was used to keep these flat and the last 300 ns of each simulation used. All standard errors are all <0.001 nm2.

Lipid tail order parameters (S) were calculated as described in Ingólfsson et al.[30] That is, the angles (θ) of each bond from the linker to the lipid tail with respect to the Z-dimension (an approximation of the bilayer normal) were evaluated and calculated. Averages of all bonds are reported as well as S between tail beads 2 and 3. All values are from 2 μs averaged at the end of each simulation, and when averaging between lipids, the averages are weighted based on the number of lipids.

The lipid lateral diffusion coefficients (D) were calculated using the GROMACS tool g_msd from the membrane plane mean square displacement (MSD) of the ROH, GL1, or AM1 beads for the cholesterol, glycerol, or ceramide lipids, respectively. The last 10 μs of each simulation was used for the analysis and straight line least-squares fitted to the MSD curve, excluding the first and last 10%, to determine D. Diffusion coefficients are reported as is and not scaled due to the faster effective dynamics of the CG force field[55] nor due to finite system size artifacts.[56]

Cholesterol flip-flop is defined and measured as described in Ingólfsson et al.,[30] where a flip-flop is counted when a lipid moves from one leaflet to another. Cholesterol flip-flop was calculated for the last 10 μs of each simulation for each lipid class using a cutoff length for what is considered within a leaflet to 1.1 nm to removing spurious flip-flop events. Reported errors are standard errors of the mean estimated by splitting the last 10 μs of the simulations in three equally sized blocks and analyzed separately.

Nonideal lipid mixing was evaluated by counting lipid neighbors for each lipid species as described in Ingólfsson et al.[30] For each lipid type, all lipids within 1.5 nm radius in the bilayer plain were counted. The last 10 μs of each simulation was averaged and, for each lipid type, the relative enrichment/depletion reported as compared to a random mixing of each mixture.

Results and Discussion

Bilayer Compositions

To generate an initial mixture for the reduction step from ∼60 to ∼15–20 lipid species, a set of guidelines are followed. The different types of headgroups and their relative populations should be maintained. Within the different headgroups, lipid species with similar tails are merged into a single representative species. These merged lipid species for the different headgroups should recapture the average number of unsaturations per tail for that headgroup. For example, the inner leaflet of A-63 contains ∼1600 phosphatidylcholine (PC) lipids (in seven different species), of which ∼500 are POPC, ∼800 are PIPC (C16:0/18:2), and the rest are mainly lipids with one polyunsaturated tail. For the A-18 composition, the seven PC species are down-selected to just ∼500 POPC and ∼1100 PIPC, a ratio that as much as possible retains average unsaturations per tail for PC lipids, as well as the distributions of unsaturations per tail. This process is repeated for all lipid classes with the goal of producing an initial mixture that reduces the number of lipid species from ∼60 to ∼15–20. The global tail unsaturation properties (average unsaturations per tail, distribution of number of unsaturations) for these mixtures are calculated. To optimize the agreement with the global tail unsaturation properties of the ∼60 lipid composition, the type and percentage of the lipid species in the initial 15–20 lipid mixture are then collectively adjusted.

The procedure is repeated to further down-select the compositions as far as possible while still being able to retain the global headgroup and lipid tail distributions. Our iterations reveal that a mixture of eight lipids appears to be an approximate threshold for making the compositions “reduced”. Attempts to generate mixtures of fewer lipid species resulted in compositions that failed to satisfactorily recapture even the basic head/tail distributions of the ∼60 lipid bilayer; the compositional flexibility is just not present.

The final distributions of the lipid headgroup classes and lipid tail unsaturations for the six different membrane compositions are displayed in Figure (full details of the exact compositions are included in the Supporting Information). A challenging task when attempting to use fewer types of lipids to reproduce the distributions seen in the most complex (∼60 lipid species) membranes is the reduced specificity in the lipid types available. One could pick lipid species that precisely reproduce the headgroup distributions but are unable to recapture the spectrum of lipid tail unsaturations. This leads to an iterative process of adjusting the composition of the ∼16 and 8 lipid mixtures to find a compromise and balance between accurately recapitulating both the headgroup and lipid tail distributions.

Figure 1

Lipid distributions of the different lipid compositions. Pie charts show the distribution of different headgroup classes between the inner and outer leaflets. The level of tail saturation for each leaflet is also illustrated.

When reducing the lipid species from ∼60 to ∼16, the proportions of the major headgroup classes can be maintained almost exactly, with only relatively minor changes to the number of lipid tail unsaturations. Indeed, the pattern of lipid tail unsaturations remains clearly distinct between the A-18 and B-16 compositions. However, when the number of lipids is reduced further to just eight in the A-8 and B-8 compositions, concessions have to be made. For instance, the “other” category of lipid headgroup (such as PI and PA lipids and ceramides) completely disappears, as there simply are not enough different lipid types. Furthermore, not even all of the main headgroup classes can be maintained and selection of which to keep should be based on what is the target application for that mixture; e.g., the outer leaflet glycolipid component in the A-63 and A-16 compositions (which comprises only ∼5% of the lipids in the leaflet) can no longer be accommodated in the A-8 mixture.

However, replicating the average distribution of categories of lipids does not guarantee replication of how the lipids will behave together within the membrane. For instance, a membrane of 80% POPC and 20% cholesterol would have the same average unsaturations per lipid tail (0.5 unsaturations per tail) as a membrane of 40% DOPC, 40% DPPC, and 20% cholesterol, but their behavior would be drastically different. In order to compare the behavior of the systems, global membrane and leaflet analysis is required.

General Bilayer Properties

Multiple general bilayer properties are calculated and compared for all six of the CG simulation systems (Figure ). The two most significant features that we want to determine are (a) if these properties are maintained as much as possible compared to the “gold standard”, highly complex ∼60 lipid species simulation and (b) if the properties begin to deviate, there is still a distinct difference between the membranes of different tissue type.

Figure 2

Global membrane properties. For both the inner and outer leaflets of the six membrane systems, the average number of unsaturations per lipid tail (A), cholesterol fraction in each leaflet (B), average area per lipid (C), average lipid order parameter at position 3 of the lipid tails (D), and average lipid diffusion rates (E) are calculated. Where applicable, the standard error is shown. All of the data plotted is also reported in Tables S2 and S3.

Several properties are found to be extremely consistent for the different simulations of each tissue type. For both leaflets, both tissue types, and both the ∼16 and 8 lipid complexity, the cholesterol fraction (Figure B) and average area per lipid (Figure C) display <3% deviation when compared to the ∼60 lipid simulations. This also importantly means that these properties remain distinct between the “Average” and “Brain” simulations even when the lipid complexity is reduced to ∼16 or 8 lipids. In fact, almost all of the properties measured remain different between the average and brain simulation types. The only exception to this is the average tail order parameter for the inner leaflet of the B-16/8 and A-18/8 systems. Here, the B-16/8 order parameters are more consistent with the A-63 system, while the A-18/8 order parameters are more consistent with the B-58 system. Despite this overlap of order parameters, the equivalent areas per lipids for these leaflets are almost identical to their comparative ∼60 lipid tissue simulations, and the order parameters for the outer leaflets of these simulations remain dissimilar from each other. One may attribute this deviation in order parameter to the marked drop in average unsaturations per tail of the chosen B-16 inner leaflet composition (Figure and Figure A, pale red line). However, the order parameter has contributions from many other properties such as the cholesterol fraction and the lateral spatial distribution of the lipids. This again highlights that balance that must be achieved when attempting to reproduce the overall behavior of a bilayer, rather than specifically overfitting to a single property.

Additional properties are calculated for the membrane as a whole, rather than individual leaflets (Table ). As cholesterol can flip-flop between the two leaflets during the simulation time scale, a transition rate can be calculated. The cholesterol flip-flop rate is shown to be consistently higher for the A-63/18/6 membranes (∼(6–8) × 106 s–1) than for the B-58/16/8 system (∼(4–5) × 106 s–1). With A-8 somewhat surprisingly lower than A-63 and A-18. There is a clear difference between the cholesterol flip-flop behavior of the different tissue types, regardless of the complexity of the membrane. The bilayer thickness, however, is not as clear, with the Brain mixtures tending to be marginally thicker. Relatively large fluctuations within the thickness of the membrane (as well as potential in precisely calculating the thickness) mean that any trends are challenging to identify, with all systems within (or close to) standard deviation.

Table 2

Membrane Properties

system	A-63	A-18	A-8	B-58	B-16	B-8
bilayer thickness (nm) ± SD	4.109 ± 0.064	4.105 ± 0.026	3.942 ± 0.019	4.137 ± 0.143	4.182 ± 0.047	4.131 ± 0.026
CHOL flip-flop (106 s–1)	7.290	7.700	5.910	4.820	4.373	4.374

Lipid Neighbor Analysis

Having quantified that the tissue simulations of different complexity have very comparable global behavior, we sought to investigate some of the more local properties. The nonideal mixing of lipids with different headgroups is classified by calculating the local enrichment or depletion of different lipid species (Figure ). The local lipid analysis reveals a complex pattern of behavior, for which parsing out information can be difficult. There are some general trends seen within the Martini force field that are consistent between different tissue mixtures and are maintained throughout the differing levels of complexity; glycolipids (where present) and PIPs both self-associate, while PIPs also have decreased interactions with the other anionic lipids (PS), and sphingomyelins (SM) show decreased association around PE headgroup lipids. Mutual coordination of counterions mediates the aggregation of PIPs,[57] while electrostatic repulsion likely accounts for the decreased association between the different anionic species. The decreased interaction of SM around PE headgroup lipids appears to be due to the distinct differences in the tail saturations of these lipids. The SM lipids consist mainly of saturated hydrocarbon tails, such as palmitoyl (C16:0), arachidoyl (C20:0), or lignoceroyl (C24:0). In contrast, the PE lipid tails have a high proportion of polyunsaturated tails, such as arachidonoyl (C20:4) and docosahexaenoyl (C22:6). Conversely, due to the more saturated nature of their tails, SM lipids consistently show enrichment of cholesterol around them (particularly in the inner leaflets). There are also some subtle differences in enrichment patterns between the tissue types; there is a stronger self-interaction between PE headgroup lipids in the Brain mixtures than the Average ones (specifically in the outer leaflet), and Brain glycolipids display significantly reduced interactions with PE headgroup lipids. This difference may be due to the Brain glycolipids containing a substantial fraction of cerebrosides, whereas all of the Average glycolipids are gangliosides. The glycolipid behavior is quite consistent within the different levels of Brain tissue complexity, as are most of the predominant trends. However, while the trends remain similar, some of the local enrichments/depletions become enhanced in the B-8 system. This is likely due to the reduced spectrum of lipids available, particularly with respect to the lipid tails. Again, using the PE headgroups as an example, the B-58 mixture has a wide array of six different PE tails, with various combinations of 1, 3, 4, 5, and 6 total unsaturations (a mean value of 4.29 unsaturations per lipid). While the B-16 mixture can still span most of this range with 1, 4, and 5 unsaturations (and a mean of 4.11 unsaturations per lipid), B-8 only has a single PE lipid with 4 total unsaturations to best match the mean unsaturations per lipid. A further example is the reported enhancement in SM–SM enrichment seen in the B-8 mixture compared to the B-58 and B-16 mixtures. The SM lipids in B-58 are actually a mixture of DPSM (d18:1/18:0), POSM (d18:1/18:1), PNSM (d18:1/24:1), and PBSM (d18:1/C22:0); B-16 is DPSM and PNSM. However, B-8 contains only DPSM. Given that the SM lipids of B-58 and B-16 contain a mixture of completely saturated and monounsaturated lipid species, it is perhaps unsurprising that they exhibit lower SM–SM association than the B-8 SM category, which is entirely DPSM.

Figure 3

Local enrichment and depletion of different lipid classes. The normalized number of lipids of a specific class (columns) within 1.5 nm of other lipid classes (rows). The data is grouped by lipid headgroup type. Enrichment is indicated with blue colors, while depletion is highlighted in red.

Domain Formation and Behavior

A further important property of a membrane that is simultaneously a combination of both global and local behavior is the partial demixing of lipids into laterally heterogeneous regions, or domains. We have previously demonstrated that the tissue-specific membrane simulations (A-63 and B-58) have notably different behaviors in terms of the size, number, and composition of their cholesterol-enriched and cholesterol-depleted domains.[31] The domain sizes are affected by bilayer undulations[31] and the domain register between leaflets is affected by cholesterol flip-flip.[64] The very generalized difference between the tissue types is that the A-63 system has larger domains, while the B-58 system has domains that are smaller, more numerous, and generally contain higher cholesterol density.

Due to cholesterol being the most consistent single factor between all levels of membrane complexity, it is again used as the metric for defining domains. The local cholesterol distribution function in the leaflet is calculated (as described in the Methods) and used to define the regions with either a “depleted” or “enriched” cholesterol density (i.e., domains). In order to define these domains, we determined a depleted and enriched threshold of cholesterol density (cholesterol distribution function) using the same methodology as that in Ingólfsson et al.,[31] where the distribution of the number of domains with respect to the density threshold was noted (Figure S1). If the choice of the enriched threshold is too high, then this stricter criterion will result in some enriched regions not qualifying as domains. If the choice is too low, then several distinct regions of enrichment may merge into a single domain. The reverse argument is also true for the choice of the depleted cholesterol threshold (choose too high and regions will merge, choose too low and domains will be missed). As such, the thresholds are chosen as the peak values of the distributions of domain counts (Figure S1). The calculated thresholds (per leaflet) to determine the cholesterol-enriched and cholesterol-depleted domains for each simulation system are shown in Figure , with a visual example of the domains illustrated in Figure S2.

Figure 4

Cholesterol-enriched and -depleted domain thresholds. The thresholds for the cholesterol-enriched and cholesterol-depleted domains are shown for the six simulation systems. The orange/gray boundary indicates the depleted-cholesterol threshold, and the gray/purple boundary indicates the enriched-cholesterol threshold. The bars above the dashed line are the outer leaflet thresholds, while the ones below the dashed line are the inner leaflet thresholds.

As to be expected, due to their inherently higher cholesterol content, the three Brain systems have thresholds that are higher than those of the Average systems. Furthermore, as the inner leaflets have a greater fraction of the cholesterol, they also have slightly higher thresholds. The threshold levels themselves are extremely consistent within the Brain systems. While more variability is seen within the Average systems, they still exhibit behavior much closer A-63 than any of the Brain compositions.

Using these system-specific thresholds, each simulation is analyzed to produce the distribution of the number of domains (per μm2) with respect to the mean cholesterol distribution within that domain (Figure ). There are two noticeable trends within the distribution of the number of domains. First, all domains for the Brain systems have higher mean cholesterol density than the equivalent Average systems (again, as to be expected, given the increased cholesterol fraction in the Brain compositions). Second, there is strong internal uniformity within each separate category (tissue type, enriched and depleted domains, inner and outer leaflets), emphasizing a consistency in the domain behavior across tissue simulations of differing lipid complexity. Additionally, there are steady trends in the mean domain size (in terms of area; nm2), whereby the domains for the Average systems are larger than the domains for the Brain systems (Figure S3). As with the domain density, there is also strong internal consistency in domain size between the tissue systems of different lipid complexity (with A-8 being the only outlier, though this may be related to finite size effects).

Figure 5

Distribution of domain cholesterol densities. The distribution of depleted (left) and enriched (right) domains for the various Brain (top) and Average (bottom) systems. The solid lines represent the outer leaflet, and the dashed lines represent the inner leaflets. Due to the different system sizes, the number of domains has been normalized per μm2.

There is one further interesting feature of note that was quite consistent between all levels of complexity; there appears to be overlap (or approaching overlap) between the cholesterol density in the enriched domains of the Average compositions and the depleted domains of the Brain compositions (Figure S4). This leads to a potentially fascinating hypothesis, that, despite their differing lipid mixtures, there may be subregions within the Average and Brain membranes that actually have remarkably similar compositions (as governed by the cholesterol density). This overlap of lipid compositions (and thus local membrane behavior) would allow the same membrane protein to exist within both tissue bilayers and experience similar lipid environments/properties.

All-Atom Simulations

The previous analysis establishes that most of the global bilayer properties, local lipid enrichments/depletions, and domains of differing cholesterol density are quite well maintained as the variety of different lipid species is reduced from ∼60 to 8. While not every feature is preserved from their parent A-63 and B-58 mixtures, the A-8 and B-8 systems are significantly distinct in their behavior and retain an individual character that identifies them as one tissue type or the other. Given that this tissue-specific behavior can be observed using eight lipid species at the CG level, converting (“backmapping”[45]) the systems to AA representations acts as a proof-of-principle test and verifies at a high level that these systems are stable and retain the broad features of the membrane.

The ∼6000 lipids of the A-8 and B-8 systems are still quite large for AA simulations. As such, the equilibrated coordinates from A-8500 and B-8500 (the same compositions as A-8 and B-8 but only ∼500 lipids in size) are used to generate the starting configurations for the equivalent AA simulations, A-8AA and B-8AA (Figure ).

Figure 6

Coarse-grained and all-atom snapshots. The equilibrated final coordinates from the CG simulations and the equivalent generated starting configurations for the AA systems. The headgroup and tail colors are as in the pie charts in Figure .

Density profiles comparing the peaks of various components of the membranes (such as lipid headgroups, linker regions, and tails) for both the A-8500/A-8AA and B-8500/B-8AA systems show equivalent patterns for both sets of simulations (Figure S5), as well as comparable differences between the Average and Brain systems. Consistent with our earlier work,[31] both B-8 systems show slight broadening of the lipid headgroup peaks compared to the A-8 systems, while cholesterol density is also increased within the middle of the B-8 bilayers. Notably, the AA systems are slightly thicker (∼2 Å) than their CG equivalent simulations. This difference is consistent with the specific choice of representative lipid tails for the AA systems (see Methods and Table S1) compared to the CG building block approach where each tail represents a range of AA lengths. For instance, many of the palmitoyl tails (C16:0–C18:0 in Martini) in B-8500 are converted to the more biologically common and slightly longer[58] stearoyl tails (C18:0) in B-8AA.

The A-8500 and B-8500 simulations were extended a further 25 ns from the frame that was used to generate the A-8AA and B-8AA systems. Due to the approximate 4-fold increase in CG dynamics,[32] these 25 ns of CG simulation were evaluated against the first 100 ns of A-8AA and B-8AA where the lateral 2D density of cholesterol, saturated tail regions, and unsaturated tail regions were compared (Figure S6). The enrichment/depletion patterns of these density plots indicate a spatial consistency of the lateral distributions between both the AA systems and CG simulations they were generated from. Thus, these relatively short, detailed atomistic simulations act as a proof-of-principle to demonstrate the analogous representation of their eight-component CG counterparts and, by extension, comparable to their highly complex ∼60-component CG systems.

Conclusions

This work has shown that many composition-dependent, tissue-specific properties of CG bilayers, both general and specific, can be retained when the complexity of that composition is reduced. Furthermore, we can approximate the behavior of the CG membrane with reasonably reduced number of components to allow for detailed AA simulations. However, not all properties are retained, especially with the A-8 and B-8 systems that only contain eight components. There are trade-offs and concessions that have to be made when the dynamic plasticity of the membrane elements is reduced to just eight different lipid species (and really just seven, if cholesterol is not counted).

Thus, if wanting to use a reduced eight-component system, one should take into account what is the focus/purpose of those simulations. Are there certain membrane properties of interest? Or is the membrane being used to study a specific protein? The composition of the system may have to be adjusted depending on the type of scientific question being asked. Unfortunately, there does not appear to be a formula or algorithm to determine which mixture will be most appropriate. In the A-8 mixture, the absence of glycolipids is conceded, as they only account for ∼5% of the headgroups in the A-63 outer leaflet. However, if the goal of the study is to investigate the behavior of a specific membrane protein that may interact with glycolipids, then evidently the composition would have to be reassessed. Indeed, PIPs are included in system A-8. While only accounting for ∼2% of the inner leaflet, PIPs are known to be key lipids for interacting with certain proteins, and system A-8 was designed with that in mind.[59] Conversely, by limiting ourselves to just eight lipid types and deciding to include a glycolipid species in B-8 (as there are >10% glycolipids in the B-58 outer leaflet), it means compromising elsewhere by only having a single type of PE lipid. What is gained in terms of richness of headgroup diversity is relinquished in range of lipid tails.

There have been “complex” (asymmetric bilayers with greater than five different lipid species) mammalian membrane model systems that have been utilized to study several aspects of membrane and membrane protein behavior to great effect.[17,60−63] However, through direct comparison with systems of varying complexity, we are able to validate our two different eight-component bilayers that can act as reliable mimetics for our biologically complex Average or Brain membranes. Using systems of reduced complexity and smaller size significantly lessens both the computational costs and the complexity of analysis. These systems also allow for equivalent study using all-atom simulations and even artificial membrane experiments.