Modeling lipid raft domains containing a mono-unsaturated phosphatidylethanolamine species

Abstract

Several membrane proteins are preferentially partitioned in lipid microdomains called rafts. The hypothesis of an intimate relationship between proteins and their specific raft environment is nowadays widely accepted. Indeed, the raft–protein cross-talk would influence protein activity and trafficking either by specific lipid–protein interactions or changes in physico-chemical properties of the bilayer. Although lipid rafts used to be simply considered membrane patches enriched in sphingolipids, cholesterol, and saturated phosphocholine derivatives, the optimization of extraction procedures and recent lipidomic analyses challenged this established concept, highlighting a significant presence of phosphatidylethanolamine species. Relying on this evidence, we devised a generic coarse-grained raft-like model containing di-stearoyl phosphatidylcholine, cholesterol and palmitoyl-oleoyl phosphatidylethanolamine species. The model was validated against available experimental data by studying the lipid mixture at different molar ratios through extended molecular dynamics simulations. The agreement of structural and dynamical properties with those of a liquid-ordered crystalline phase suggests that our model can represent a reliable lipid environment especially suited for computational studies aimed at unraveling raft–protein functional interactions.

---

1. Introduction

Under physiological conditions, the components of biological membranes can group according to their physico-chemical properties, achieving a non-homogeneous lateral distribution and forming distinct domains. These domains exist in a more structured and ordered fluid phase (liquid-ordered, Lo) floating and diffusing in the surrounding liquid-disordered (Ld) phase.^1,2 In living cells, these regions are called lipid rafts and are conventionally identified as “small, (10–200 nm) heterogeneous, highly dynamic, sterol- and sphingolipid-enriched domains that compartmentalize cellular processes”.³ An ever increasing number of membrane proteins such as tyrosine kinases,⁴ GPCRs and neurotransmitter transporters have been reported as being preferentially partitioned in lipid rafts.^5–7 Although lipid membranes have traditionally been regarded as passive solvents, it is nowadays recognized that lipids have biological functions and play a key role in several biological events. Signal transduction, membrane trafficking, and protein conformational equilibrium, are examples of processes that can be modulated by the lipid composition of the membrane.² In turn, by preferentially interacting with raft-like lipids, proteins can facilitate the assembly of the raft itself.⁸

The nature of lipid rafts was originally deduced by the composition of the insoluble fractions obtained by treating biological membranes with the non-ionic detergent Triton X-100.¹ Resistance to such detergent is observed in bilayers showing a tight lipid packing, like domains in the Lo phase. The resulting insoluble fractions, called detergent-resistant membranes (DRMs),⁹ were found to be enriched in sphingolipids (SLs), phosphatidylcholines (PCs), and cholesterol (Chol) and depleted in glycerophospholipids (GPLs) such as phosphatidylethanolamines (PEs) or phosphatidylserines (PSs). While saturated SLs and PCs reside preferentially in the extracellular side of the membrane, where they act as stable physical barrier, PEs and PSs are typically located in the cytoplasmic leaflet, and are characterized by inability to form bilayers by themselves.¹⁰ These evidences led to the conclusion that lipid rafts must be selectively enriched in Chol and SLs over GPLs and had to exist only in the outer leaflet of membranes.¹ This body of evidence, however, was severely challenged by comparative experiments making use of milder detergents which demonstrated that lipid rafts could also be found in the cytoplasmic side, even though with different lipid compositions and phase properties.^10–13 In particular, DRMs obtained by detergents belonging to the Brij series retained the physiological asymmetry in phospholipid distribution of membranes, and contained fully functional membrane proteins.^{4,6,7,14–16} These so-called “atypical” DRMs (to be readily distinguished by “traditional” DRMs) displayed a lower level of Chol and SLs and, consistently, higher concentrations of GPLs with higher fractions of unsaturated species. Recently, the discrepancies between traditional and atypical DRMs have been addressed by Morris and co-workers, who showed that by optimizing experimental conditions, a comparable enrichment in GPLs, SLs, and Chol could be obtained.^17,18 In spite of this, the fact that DRMs encompass multiple compositions is an indirect evidence that in cells lipid rafts must be a collection of highly heterogeneous domains, whose physico-chemical properties are finely regulated by changes in relative abundance of lipid components within the domain itself.¹⁹

From a computational standpoint, both atomistic (AA) and coarse-grained (CG) Molecular Dynamics (MD) simulations have been extensively used to explore physics and biology of lipid rafts. Although the former provides a highly detailed picture of interatomic interactions, CG-MD, exploiting a particle-based representation, allows one to observe events naturally occurring on time- and length-scales that are typically out of reach for standard AA-MD, such as lipid diffusion and phase-separation.²⁰ Bilayers including typical raft-like components, such as saturated PCs (or SLs), poly-unsaturated PCs, and Chol, have been extensively employed within a CG framework for their ability to phase separate in Lo and Ld domains, the former being enriched in Chol and saturated species, and the latter mainly composed by the poly-unsaturated ones.^21,22 Moreover, it has been shown that compositional and phase asymmetry could be reproduced with these models, simulating symmetric bilayers in an homogeneous Lo state, as well as asymmetric patches consisting in a Lo leaflet coupled to a Ld one.²³ Even though the properties of Lo phases have been largely investigated for lipid species belonging to the outer leaflet, no CG-MD studies of inner leaflet lipids in a pure raft-like assembly have been reported so far. On the contrary, concerning AA-MD simulations, Bhide and co-workers have modelled a completely asymmetric raft-like system containing mixtures of SM/Chol in the outer leaflet and an unsaturated PS derivative/Chol mixture to mimic the cytoplasmic leaflet.²⁴ More recently, a series of AA-MD simulations has been performed to assess the structure and dynamics of symmetric raft-like mixtures surrounded by a completely asymmetric Ld bulk phase.²⁵ While this study took into account the contribution of unsaturated species in Lo domains, inner leaflet lipids were represented only in the Ld phase. Given the relatively limited time-scales of the simulations performed and the variability of lipids involved in rafts, there is still a compelling need to investigate these assemblies by validating reliable raft-like lipid mixture models.

In line with latest lipidomic evidences, we devised a generic CG model for an atypical raft bilayer characterized by the presence of a lipid species representative of the inner leaflet of biological membranes. Specifically, we assembled a three component bilayer containing palmitoyl-oleoyl phosphatidylethanolamine (POPE), di-stearoyl phosphatidylcholine (DSPC) and Chol at three different molar ratios with the MARTINI force field.²⁶ To the best of our knowledge, this is the first ever reported CG investigation of a pure raft-like model containing mono-unsaturated PE species. We show that through a careful selection of lipid components, in each case, a homogeneous Lo phase was observed and preserved throughout 10 μs of CG-MD simulation. The ability of such lipid mixture to display structural and dynamical properties in agreement with those reported for biologically relevant Lo phases makes it an advanced model able to reliably capture the salient features of lipid rafts in computer simulations. Moreover, being purposely designed to be of general applicability, such a model is expected to be useful to study raft–protein interactions whenever specific lipidomic data for the system under investigation is missing. Indeed, provided that conformational changes in membrane proteins can only be reproduced by explicitly simulating the specific lipid framework in which the protein is embedded, this work may represent a step toward a complete customization of the lipid environment for the simulation of membrane-protein systems.

2. Computational methods

Lipid selection and model system

Even though the composition of biologically relevant lipid rafts is highly heterogeneous, the lipid mixture employed in our model was specifically chosen so as to exhibit a Lo phase in standard temperature and pressure conditions. This requirement was obtained by identifying some common features of lipid rafts that were exploited to design a generic model of an atypical raft domain.¹⁹ In particular, (i) the enrichment in outer-leaflet lipids containing long and saturated acyl chains, (ii) a moderate enrichment in Chol, and (iii) the presence of inner-leaflet mono-unsaturated species were explicitly taken into account. Thus, relying on experimental and theoretical data, an atypical raft model consisting of DSPC, POPE, and Chol was generated. We note that, differently from several previously reported CG models, we were mostly interested in investigating and reproducing stationary properties of lipid-rafts, rather than directly observing the Lo/Ld phase separation.²²

Phospholipids or SLs with high melting temperatures can tightly pack in the presence of Chol to form Lo domains and phase-separate from a large variety of unsaturated PCs with lower transition temperatures. In particular, DSPC displayed the ability to form Lo phases comparably to SL species such as sphingomyelin (SM) when mixed with di-unsaturated PCs, PEs, and Chol at the temperature of 298 K, well below its melting point (T_m,exp = 328 K).²⁷ For this reason, fully saturated PC derivatives have systematically been employed as a substitute of SLs to reproduce the formation of raft-like phases in model membranes.^28,29 We note that SM and DSPC share the same head-group and the same stearic acyl chains (C18:0). The latter, together with palmitoyl chains (C16:0), were shown to be well represented in raft-related components (GPLs and SLs derivatives) and to have a key role in promoting and stabilizing raft formation via favourable interactions with Chol.^1,2,19,27 Given that PCs are the major components of biological membranes, accounting for approximately 50% of total lipids and are equally distributed between leaflets,^12,30 we decided to use DSPC instead of SM as a representative of the high melting lipid in our raft-like model. Similar considerations led us to identify POPE as a candidate to represent GPL species in atypical raft domains. Although PE-Chol interactions are in general considered as relatively unfavourable, recent studies revealed that some mono-unsaturated PE species do not phase-separate from Chol-enriched domains,^31,32 and they can also be able to accommodate up to 52 mol% of Chol.³³ These considerations made us confident in using POPE as second raft-like lipid for our model. To take into account a moderate enrichment in Chol as found in atypical rafts, a maximum sterol concentration of 25% was considered.

The parameters employed to model the selected lipids were taken from the MARTINI force-field v2.0.^26,34 MARTINI uses a 4 to 1 mapping scheme whereby on average four atoms are represented by a single interaction bead particle, except for ring-like groups that are mapped with a 3 to 1 strategy. In details, the DSPC head-group consisted of two hydrophilic charged groups: the choline (type Q₀) and the phosphate group (type Q_a). The positive ethanolamine in POPE was modelled through a Q_d bead particle. In both lipids the glycerol group was represented with two sites of intermediate hydrophilicity (N_a). Each of the stearoyl tails was mapped into 5 apolar beads (C1), whereas palmitoyl and oleoyl chains were modelled with 4 beads (C1), and 5 beads (C1 and C3), respectively, where the C3 bead-type was used to take explicitly into account the polarizable nature of the double bond in oleoyl tails. Chol was modelled through 8 bead particles mapped on a 3 : 1 basis in order to represent fused rings.

Six of these beads (ROH, R1–R5) represented the sterol body and the remaining (C1, C2) the short tail. A proper combination of bond, angles, and dihedral constraints was used to maintain the average planarity of the fused rings of Chol (see Fig. 1).

Fig. 1 Atomistic (right) and coarse-grained (left) representation of the lipids used in the raft-like membrane model. For clarity, only polar hydrogens are shown in the atomistic representation.

Setup of the systems and molecular dynamics simulations

To investigate the condensing effect of the sterol on lipid properties, three raft-like membranes with an increasing amount of Chol were built. These final systems were obtained by as many homogeneous patches generated with the MARTINI self-assembly protocol³⁵ and contained 218, 204 and 192 DSPC molecules (denoted as patch 1, 2 and 3, respectively). In details, for each patch, DSPC molecules together with 6 CG water beads per lipid (corresponding to 24 effective water molecules per lipid) were randomly mixed in a rectangular box, minimized, and equilibrated in the NPT ensemble under isotropic pressure conditions at 300 K and 1 bar. MD simulations were performed using periodic boundary conditions with a neighbour list cut-off of 1.4 nm and a time-step of 20 fs. Non-bonded interactions were treated accordingly to the MARTINI model standard protocol.²⁶ Temperature and pressure were controlled using the Berendsen weak coupling algorithm.³⁶ A compressibility of 3 × 10⁻⁴ bar⁻¹ and a relaxation constant of 5 ps were used. All simulations were performed with Gromacs 4.6.³⁷ Because of the speed-up factor due to the CG representation, effective simulation times are approximately obtained multiplying by a factor of 4 the amount of actual sampling.²⁶ For clarity, throughout the paper results are reported over the actual simulation time, unless explicitly stated. The three patches were simulated for 200 ns. While patch 1 and 2 returned defect-free DSPC bilayers at 300 K, in order to avoid phase transitions, patch 3 was simulated at 330 K (above macroscopic melting temperature, T_m,exp).²⁷ At the end of the self-assembly procedure, 90, 76, and 64 DSPC molecules were replaced with POPE molecules in patch 1, 2, and 3, respectively. In this way, each patch contained the same number of DSPC molecules, which mimics the major lipid components (SLs and PCs) in our raft-like membrane models. The stress introduced by molecule replacement was reduced by an additional minimization. Then, the three patches were further equilibrated for 200 ns switching from isotropic to semi-isotropic pressure coupling. Finally, 38, 52, and 64 Chol molecules were introduced in patch 1, 2, and 3, respectively, and the systems were thus energy minimized and equilibrated once again for 200 ns. The molar ratio of Chol was increased up to the limiting solubility of this molecule in POPE, corresponding to 30, 40, and 50 mol% relatively to the already present amount of POPE in each patch, and hence to a total of 15, 20, and 25% of Chol. These ratios were chosen so as to optimize the DSPC/Chol mixing,³⁸ and also because of the availability of experimental data reporting the condensing effect of Chol at these concentrations.³¹ The final systems, obtained by replicating these smaller patches along axes x and y and hereafter named as Chol_15%, Chol_20%, and Chol_25% (see Table 1), consisted of a total number of 1024 lipids, and were then simulated for 200 ns before acquiring averages. All production runs were performed at 300 K, which is close to DSPC transition temperature observed in previous simulations (T_m,calc = 295–299 K).^39,40

Table 1 Lipid composition of the simulated systems (molar ratios)

	Chol15%	Chol20%	Chol25%
DSPC : POPE : Chol	0.50 : 0.35 : 0.15	0.50 : 0.30 : 0.20	0.50 : 0.25 : 0.25
POPE : Chol	0.70 : 030	0.60 : 0.40	0.50 : 0.50

---

Since DSPC was also selected to mimic SLs, this choice was instrumental to capture physiological aspects that can be important to build a well-behaved raft-like model. Indeed, in biological membranes, SLs are the species with higher transition temperatures (above 310 K), and they are very concentrated in many important tissues where they show a liquid crystalline phase.⁴¹ It has been suggested that this property can be, at least in part, responsible for the recruitment of raft components within a Lo environment.² The initial structure of the raft-like membrane systems are reported in Fig. 2. Each system was simulated for 10 μs.

Fig. 2 Top and side view of the membrane models at the end of the equilibration. DSPC, POPE, and Chol are represented in blue, red, and green, respectively.

Analysis of trajectories

MD trajectories were analysed to extract structural properties and long-time diffusion profiles for single species. Individual areas per lipid were obtained via a grid-based analysis, averaging in time values extracted every 200 ps for single leaflet with the GridMAT-MD software.⁴² The bilayer thickness was calculated in a similar way with the same software.

The structural order of the system was assessed by extracting the deuterium order parameter (S_n) from trajectories via the do-order-multi.py script available with MARTINI.²⁶ The order parameter is defined as follows:

S_n = 1/2(3〈cos² θ_CH〉) − 1 (1)

where θ_CH is the angle between a CH bond and the bilayer normal. S_n equals to 1 means a perfect alignment with the bilayer normal, values approaching −0.5 result in anti-alignment, and S_CD = 0 indicates a random orientation. Within a CG description, the vectors used to calculate the angle lied on the direction defined by the bond between two adjacent beads. The averages were collected on time windows of 1 ns each.

Diffusion in the xy plane was obtained by fitting a straight line to the time-dependent Mean Squared Deviation (MSD) of atomic positions. The 2D Einstein's relation establishes a direct connection between two-dimensional diffusion constant and the limiting slope of the MSD profile:

here, r_i(t) is the position of lipid i at time t, and r_i(0) is the position at which the measurement starts. The lateral diffusion constants were calculated from the average lateral MSD of a representative bead in the head group region for each lipid.⁴³ The phosphate bead group (PO4) was chosen for DSPC and POPE, while the hydroxyl bead group (ROH) was employed for Chol. Single-leaflet calculations were performed and then averaged to obtain a unique MSD profile for each species in each system. For each simulation, the diffusion coefficient was calculated averaging the unique MSD profile over 10 evenly distributed chunks along the trajectory (i.e. average over 1 μs).

The restart time for the MSD measurement was set to 50 ps. Jumps of molecules across periodic boxes were removed for the entire system, while centre-of-mass (COM) motion was removed for lipids of the same type in the considered monolayer. Because of inter-leaflet flipping, only Chol molecules remaining in the considered leaflet along the 1 μs chunk were used to calculate the MSD. Gromacs 4.6 (ref. 37) was utilized to extract MSD, as well as to calculate the density profiles. The degree of lipid mixing in the raft model was measured by calculating the preferential partitioning of DSPC and POPE in each leaflet. The fractional interactions were derived as the relative number of contacts between a pair of lipid species, normalized by the total number of contacts for all the lipids in the simulated systems:⁴⁴

Contacts were defined with respect to GL1 and GL2 beads of DSPC and POPE within a defined cut-off radius. First and second lipid neighbour analyses were carried out considering a value of 0.8 and 1.1 nm, respectively. These distance cut-offs have been used in previous works to define lipids in contact and roughly correspond to the first and the second solvation shells detectable in radial distribution functions for CG particles.^44,45 If more than one bead of the same lipid was located within the cut-off, only one contact was considered. Applying eqn (3), a fully random ideal mixture of two lipid types would result in an equal partitioning of all the components (p_x = 0.50). The analysis was performed using only the last microsecond of simulations, and repeating the calculation with a 1 ns interval. The fractional interactions were separately computed for each bilayer leaflet and results were then averaged.

The translocations of Chol molecules along the membrane normal, including complete flip-flop events, were analysed by monitoring the projection of the centre-of-mass of the sterol's body along the z-axis of the simulation box. The average equilibrium position and the inter-leaflet region were defined as layers satisfying the conditions |z| ≥ 15 Å and |z| ≤ 5 Å, respectively. To reduce the noise due to fluctuations, a translocation event was considered as accomplished when a Chol molecule, previously entered in one of the two layers, was found in the other one after a certain simulation time.

3. Results and discussion

Structural properties

In Fig. 3, the individual areas per lipid for DSPC and POPE are reported against simulation time for each simulated system, whereas the area per lipid for Chol is shown in Fig. S1 in the ESI.† As the plots show, all areas are distributed around equilibrium values indicating that systems were properly equilibrated in production runs. In systems Chol_15%, Chol_20%, and Chol_25%, average DSPC areas constantly decreased with values of 56.9 ± 0.5, 54.4 ± 0.6, and 51.9 ± 0.5 Ų, respectively. These values properly match with those obtained from AA simulations of pure DSPC bilayer performed at 300 K (ref. 39) and experiments conducted at 330 K (above DSPC phase transition).⁴⁶ Moreover, the observed decrease is in the range of a fluid phase, being quite far from values reported for gel phase (47.3 Ų).⁴⁷ While there is no significant difference between DSPC area in Chol_15% and in pure bilayers at 300 K (ref. 39) the condensing effect could be easily observed in the remaining systems, where the amount of Chol approximates the range of enrichment found in cells for atypical raft membranes. Analysing the same parameters for POPE, we found average values of 53.4 ± 0.7, 51.1 ± 0.8, and 49.0 ± 0.8 Ų for Chol_15%, Chol_20%, and Chol_25%, respectively, in line with experimental values obtained for POPE monolayers in the presence of 30, 40, and 50% of Chol at 308 K.³¹ Moreover, Chol showed an average area of about 38 Ų in all systems, displaying a very good agreement with that observed in fluid phases.⁴⁸

Fig. 3 Averaged areas per lipid for DSPC and POPE plotted as time series (left) and as probability distributions (right).

A closely related structural property of lipid bilayers is the membrane thickness. By analysing our simulations, a slight increase in thickness was found as a function of Chol concentration (Fig. S2†). In particular, we observed an average thickness of about 48 Å, which is indeed in agreement with experiments, where differences between Lo and Ld are in the order of 6–8 Å.⁴⁹ It has been previously shown that it is possible to detect such differences with CG-MD models, as an approximate thickness around 45–50 Å for Lo phases and 35–45 Å for Ld phases have been reported.^22,50

In Fig. 4, density profiles for the three raft-like membranes are reported. As it can be seen, the chemical groups show a homogeneous distribution along the bilayer normal of the membrane for all the simulated systems. Notably, despite DSPC molar ratio was kept constant among the three systems, an increasing density for DSPC acyl chains in the middle of the bilayer could be observed with increasing Chol concentrations. This indicates chains interdigitation, which has been proposed as the main factor in promoting trans-monolayer coupling.⁵¹ Indeed, domain coupling could have a functional role in the coordination of both peripheral and integral raft-related proteins activity, and is expected for atypical rafts containing inner and outer leaflet species. We can speculate that interdigitation is optimized by the Chol-induced ordering effect on the saturated DSPC tails (see below). Similarly, a non-negligible density in the same area was also found for POPE, even though, in this case, differences in density are also due to changes in lipid concentration among the three systems. Another interesting feature is the increase in Chol density in the middle of the leaflets going from system Chol_15% to Chol_25%. Moreover, from Fig. 4, it can also be noticed that the Chol density in the centre of the membrane is not zero, reflecting either the ability of this sterol to switch from one side of the membrane to the other, or a certain degree of interdigitation in the inter-leaflet region of the bilayer (or both). To better investigate this aspect, we monitored the Chol balance between leaflets (ΔN_leaf) and the number of interdigitated molecules (N_in) versus time (Fig. 5). As Fig. 5 shows, each system displayed a remarkable ability to exchange Chol molecules from one side of the bilayer to the other, while maintaining on average a balanced distribution of the sterol between leaflets. In spite of this, a small number of interdigitated Chol molecules, ranging from 0 (no interdigitation) to about 5, was also detected. As the probability distributions reported in the same Figure show, there is almost an equivalent probability to observe either 0 or 1 interdigitated Chol molecules between leaflets in each system. To further characterize the dynamical behaviour of Chol molecules, we monitored the time spent by the sterol in the equilibrium position inside the leaflets and in the inter-leaflet regions between consecutive translocations (eq → in and in → eq translocations, respectively), and we reported it as a logarithmic histogram in Fig. S3.† As it can be seen in the plots, most of the in → eq translocations occurred on a nanosecond time-scale, whereas the opposite transitions are estimated to be on average two orders of magnitude slower. From the plot, it can also be noticed that going from Chol_15% to Chol_25%, not only the total number of translocations decreased, but the relative probability associated to the slower events progressively increased, as a result of the condensing effect of the sterol. More qualitatively, such an effect can also be appreciated by comparing the time series of the sterol's body projected along the z axis for the slowest translocation events observed in each system (Fig. S4†). As a whole, we counted 925, 894, and 711 complete flip-flop events for Chol_15%, Chol_20% and Chol_25%, respectively.

Fig. 4 Chemical groups-based density profiles calculated for the simulated systems and averaged over the whole trajectory.

Fig. 5 Chol balance between leaflets and number of interdigitated Chol molecules plotted as a function of time (left) and as a probability distribution (right) for the three simulated systems. Tick marks on the probability distribution plots are each 0.05 units.

One of the most striking features of atypical raft membranes is the high degree of order in phospholipid acyl chains compared to typical Ld domains, which reflects the interaction of Chol with the lipid tails. A marked increase in order was actually found going from Chol_15% to Chol_25% system, for both saturated and unsaturated lipids (see Fig. 6). In particular, in Table 2 the order parameter averaged over the whole lipid tail for all the systems is reported. Taken as a whole, Fig. 6 and Table 2 indicate an order increase which is a function of the Chol concentration in the simulated systems. We also notice that a comparable effect was experienced by saturated stearoyl/palmitoyl chains, whereas lower ordering was found for the POPE unsaturated oleoyl chain, which is however still in line with that accepted for Lo phases. Indeed, while order parameters expected for liquid disordered phases range approximately from 0.2 to 0.4 in MARTINI simulations,^22,50 liquid ordered domains are typically characterized by values between 0.4 and 0.7.^22,34,50,52

Fig. 6 Order parameter calculated for the lipid tails of DSPC and POPE for the three simulated systems plotted against time.

Table 2 Order parameter averaged in time and over the phospholipid acyl chains. In the case of DSPC, the order parameter for the two chains is further averaged to obtain a single value

	Chol15%	Chol20%	Chol25%
DSPC stearoyl chain	0.547 ± 0.008	0.576 ± 0.009	0.614 ± 0.009
POPE palmitoyl chain	0.534 ± 0.011	0.572 ± 0.011	0.615 ± 0.012
POPE oleoyl chain	0.466 ± 0.009	0.498 ± 0.010	0.534 ± 0.011

---

Dynamical properties

Compared to Ld phases, raft-like domains display a 2 to 5 fold decrease in lateral diffusion coefficient, while keeping a remarkably fluid behaviour.⁵³ The preservation of a fluid phase is a fundamental requirement in biological membranes that has to be carefully checked in order to assess the reliability of any raft-like model. In Fig. 7, the MSD calculated for the two phospholipids and Chol is shown for the three simulated systems. The decrease in MSD observed for phospholipids and Chol going from system Chol_15%, to Chol_25%, clearly indicates a progressive reduction in lipid lateral mobility. In spite of this, all the systems maintained a fluid behaviour, as reflected by the calculated diffusion coefficient. The lateral diffusion coefficient separately calculated for each leaflet and averaged over 1 μs long chunks of the trajectory are shown in Fig. S5, S6, and S7† for DSPC, POPE and Chol, respectively. As the figures show, while the diffusion coefficient of each lipid species appreciably fluctuates along the trajectories, no major drifts are found, as a consequence of the stationary state reached by the simulated raft-like environments. Moreover, as expected, both leaflets show on average the same dynamical behaviour, reflecting their intrinsic compositional symmetry that is preserved along the trajectory.

Fig. 7 Averaged MSD for each lipid species in the three simulated systems. Averages are acquired up to 1 μs of MD simulation, and the standard deviation, regularly collected along the trajectory, is also shown as error bars.

In Table 3, the diffusion coefficients averaged along the whole trajectory as well as the standard deviations are shown. The values obtained are in line with both experiments and previously reported MD simulations.^22,50,54 Indeed, diffusion constants are in the order of 10⁻⁷ cm² s⁻¹ for phospholipids in the liquid crystalline phase, whereas the gel phase is associated to values of 10⁻⁹ cm² s⁻¹.⁴⁰ Typically, the Lo phase shows diffusion coefficients in between these ranges,⁴⁰ and the order of magnitude of 10⁻⁸ cm² s⁻¹ obtained with our calculations provides further evidence for the liquid ordered phase achieved and preserved with our raft-like models.

Table 3 Averaged diffusion coefficients calculated for each lipid species of the simulated systems. Units are ×10⁻⁸ cm² s⁻¹, and all values are opportunely scaled to take into account the CG speed-up factor

	Chol15%	Chol20%	Chol25%
DSPC	6.14 ± 0.40	4.72 ± 0.22	3.51 ± 0.26
POPE	6.41 ± 0.51	4.83 ± 0.36	3.67 ± 0.39
Chol	6.95 ± 0.96	5.07 ± 0.53	3.70 ± 0.30

---

Even though these results clearly show an informative trend in the diffusion coefficient of the considered lipid species among the simulated systems, from Table 3 it can also be noticed that for each system, DSPC, POPE, and Chol displayed a value of D in general comprised within the significant range of the other components. In other words, the dynamical properties of each lipid species appear to be more influenced by the composition of the ternary mixture rather than by the nature of the molecules themselves. Such a behavior can be partially imputed to the intrinsic limitations in accuracy of the CG representation in MARTINI. Indeed, even though this force field has been proven to quantitatively reproduce structural properties such as molecular areas, the poor representation of electrostatics and polarization effects strongly limits the possibility to observe significant differences in the diffusion coefficients due to charge contributions.⁴⁵ In spite of this, we note that dynamical properties of lipid species are, at least to a certain extent, dependent on Chol ratio, whose ordering and condensing effects are satisfactorily reproduced in MARTINI. Therefore, a Chol concentration-dependent shift of D among the systems could be expected as well.

Relying on these evidences, we speculate that a higher amount of unsaturated lipids associated with lower Chol enrichment might be responsible for the intrinsically high mobility of atypical DRMs (within the limit of Lo domains), and might contribute to phase properties of rafts in a non-negligible way.

Homogeneity of the Lo phase

All the results reported so far are in good agreement with experimentally observed properties and in general with a Chol-induced formation of a fluid raft-like environment. However, the choice of the lipid components was driven not only by the requirement to obtain a Lo domain, but also to achieve a high degree of homogeneity in the same phase. The normalized fractional interactions between DSPC and POPE reported in Fig. 8 show that there is no evidence for phase separation in our systems, as all the values are peaked around the ideal value of 0.5 for both the first and the second solvation shell. Relative differences ranging from 1 to 11% with respect to this reference value suggest a non-ideal mixing behaviour, which increases with the amount of Chol in the simulated system. As a comparison, we note that values of about 0.8 for the second solvation shell have been recently reported as a signal of weak separation in similar binary mixtures.⁴⁴ Even though a fractional interaction of 0.61 was observed in the first solvation shell for POPE–POPE in Chol_25% system, (i.e. the system where the amount of Chol is closest to that expected for atypical raft domains), this value is found to be dramatically reduced when the second shell is considered. As a matter of fact, this trend is expected to be retained in the second shell in case of a significant POPE de-mixing from DSPC.

Fig. 8 Fractional interaction matrix for DSPC and POPE in the three simulated systems. Values calculated for both the first and the second solvation shells (cut-off radius of 0.8 and 1.1 nm, respectively) are reported.

The reduction up to a 2% difference relative to the ideal mixture in the second solvation shell is indeed consistent with the presence of short-ranged local clustering, which is compatible with a single Lo phase characterized by a non-ideal mixing behaviour. It has been shown that possible mechanisms of raft stabilization in a ternary raft-like mixture containing a comparable amount of Chol (25 mol%) rely on the nature of head-groups and acyl chains of the low-melting lipids (i.e. one or two unsaturated tails).⁵⁵ Results showed that thermal stability of ordered domains increased with the head-group polarity in palmitoyl-oleoyl (PO) lipids, and POPE was found to display the higher tendency to form Lo domains. Besides, the inclusion of highly unsaturated lipids promoted phase separation from Lo saturated PC/Chol domains, whereas in PO lipids this trend was less evident. Consistently a more recent work confirmed this behaviour: in ternary mixtures containing a fixed saturated PC species, Chol and variable symmetric (two unsaturated chains) or asymmetric (one unsaturated chain) PC, phase separation was enhanced by symmetric poly-unsaturated PC. On the contrary, the presence of asymmetric mono-unsaturated PC, a unique, even if non-ideal, liquid order phase could be observed.⁵² The phase homogeneity of our systems confirmed the suitability to consider POPE as the best candidate for mixing with DSPC and Chol, both because of the nature of its asymmetric tails and due to its transition temperature (298 K), definitely higher that lower-melting PC species usually employed to facilitate phase segregation.

Indeed, it has been shown that this event is favoured in high-melting/low-melting lipid mixtures, since immiscibility is enhanced by the temperature gap between components. This is the case for highly unsaturated, low T_m PC/PE species (poor packing behaviour), whereas mono-unsaturated PE/PC species (higher T_m, and better packing properties) do not favour the process, being in an homogeneous mixing.^31,32,52

Significance of the raft-like lipid species

Here we discuss on the rationale behind the choice of DSPC as a representative of the high melting lipid species in our raft like model.

Two main theories have been reported so far to explain the formation of lipid rafts. The earlier theory focuses on the role played by acyl chains interactions,^9,28 while the other puts more emphasis on the importance of lipid head-groups.¹ According to the latter model, the presence of the amide bond in SM would be responsible for a specific interaction with the hydroxyl group of Chol which, in turn, would trigger and stabilize the raft phase. In spite of this, several studies confuted such a specificity,^56,57 even though a favorable charge pair interaction involving the choline and the polar portion of Chol was observed in computer simulations.⁵⁸ Altogether, these works support a prominent role of unspecific hydrophobic and electrostatic interactions as main driving forces in the initial stages of raft formation, rather than direct H-bonding. Indeed, van der Waals interactions among lipids and between lipids and Chol seem to be the main responsible for phase behavior in SM and PC derivatives with comparable chain properties (length and saturation).^59–61 The experimental evidence that also saturated PCs, thus lacking the aforementioned amide bond, do generate Lo phases in model membranes is an indirect proof that specific H-bond interactions are not essential to drive the process of raft formation.^28,29

In this study, a high melting PC species such as DSPC was intentionally employed as a replacement for SL derivatives to obtain a model bearing a lipid composition consistent with generic biological membranes including the inner leaflet, and at the same time showing typical properties of lipid rafts. The choice of using DSPC instead of SM was ultimately motivated by the fact that not only PCs are the major components of biological membranes, but they are also equally distributed among leaflets.^12,30 On the contrary, SLs are exclusively confined in the outer leaflet of the membrane, and therefore they are unlikely to be essential for the formation of inner leaflet rafts.⁶² Even though supporting the unspecific mechanism of rafts formation is outside the scope of our work, our results further confirm the ability of DSPC to form Lo phases when mixed with PE and Chol at room temperature.

4. Conclusions

In this work, we built and simulated an advanced raft-like model within a CG representation, and we showed that it was able to reproduce properties of liquid ordered domains in agreement with experiments and previous simulations. To the best of our knowledge, this is the first computational report on a system containing Chol and saturated PC mixed with mono-unsaturated PEs in a pure raft-like lipid framework. All raft-like components were carefully chosen so as to mimic and provide a single, homogeneous liquid ordered phase, achieved with little amount of Chol, as found in atypical DRMs. Our results support experimental evidences on raft-like mixtures containing PE species, highlighting the relevance and the compatibility of inner leaflet lipids with functional membrane domains in living cells. In light of these findings, we believe that our model mirrors the intrinsic complexity of lipid rafts in cells, and it could represent a reliable framework for future studies aimed at gaining insights into lipid rafts and lipid–protein interactions. Finally, the membrane models presented in this study are available upon request.

Acknowledgements

We thank University of Bologna, CINECA, and the Istituto Italiano di Tecnologia for financial support and computational HPC resources. We also thank the anonymous referees for their suggestions that allowed us to improve the quality of our work.

Notes and references

1. K. Simons and E. Ikonen, Nature, 1997, 387, 569–572.
2. S. Sonnino and A. Prinetti, Curr. Med. Chem., 2013, 20, 4–21.
3. L. J. Pike, J. Lipid Res., 2006, 47, 1597–1598.
4. L. J. Pike, X. Han and R. W. Gross, J. Biol. Chem., 2005, 280, 26796–26804.
5. S. K. Sadiq, R. Guixa-Gonzalez, E. Dainese, M. Pastor, G. De Fabritiis and J. Selent, Curr. Med. Chem., 2013, 20, 22–38.
6. F. Magnani, C. G. Tate, S. Wynne, C. Williams and J. Haase, J. Biol. Chem., 2004, 279, 38770–38778.
7. M. E. R. Butchbach, G. Tian, H. Guo and C.-L. G. Lin, J. Biol. Chem., 2004, 279, 34388–34396.
8. R. M. Epand, Biochim. Biophys. Acta, Biomembr., 2004, 1666, 227–238.
9. E. London and D. A. Brown, Biochim. Biophys. Acta, Biomembr., 2000, 1508, 182–195.
10. P. F. Devaux and R. Morris, Traffic, 2004, 5, 241–246.
11. L. J. Pike, J. Lipid Res., 2003, 44, 655–667.
12. V. Kiessling, C. Wan and L. K. Tamm, Biochim. Biophys. Acta, Biomembr., 2009, 1788, 64–71.
13. C. Wan, V. Kiessling and L. K. Tamm, Biochemistry, 2008, 47, 2190–2198.
14. L. J. Pike, X. Han, K.-N. Chung and R. W. Gross, Biochemistry, 2002, 41, 2075–2088.
15. E. K. Fridriksson, P. A. Shipkova, E. D. Sheets, D. Holowka, B. Baird and F. W. McLafferty, Biochemistry, 1999, 38, 8056–8063.
16. G. Radeva and F. J. Sharom, Biochem. J., 2004, 380, 219–230.
17. R. J. Morris, A. Jen and A. Warley, J. Neurochem., 2011, 116, 671–677.
18. R. J. Morris, FEBS Lett., 2010, 584, 1665–1669.
19. L. J. Pike, Biochem. J., 2004, 378, 281–292.
20. S. J. Marrink and D. P. Tieleman, Chem. Soc. Rev., 2013, 42, 6801–6822.
21. W. F. D. Bennett and D. P. Tieleman, Biochim. Biophys. Acta, Biomembr., 2013, 1828, 1765–1776.
22. H. J. Risselada and S. J. Marrink, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 17367–17372.
23. J. D. Perlmutter and J. N. Sachs, J. Am. Chem. Soc., 2011, 133, 6563–6577.
24. S. Y. Bhide, Z. Zhang and M. L. Berkowitz, Biophys. J., 2007, 92, 1284–1295.
25. F. E. Herrera and S. Pantano, J. Chem. Phys., 2012, 136, 015103.
26. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. de Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.
27. A. V. Samsonov, I. Mihalyov and F. S. Cohen, Biophys. J., 2001, 81, 1486–1500.
28. D. A. Brown and E. London, J. Biol. Chem., 2000, 275, 17221–17224.
29. D. Marsh, Biochim. Biophys. Acta, Biomembr., 2009, 1788, 2114–2123.
30. J. M. Boon and B. D. Smith, Med. Res. Rev., 2002, 22, 251–281.
31. S. R. Shaikh, M. R. Brzustowicz, N. Gustafson, W. Stillwell and S. R. Wassall, Biochemistry, 2002, 41, 10593–10602.
32. M. Grzybek, J. Kubiak, A. Łach, M. Przybyło and A. F. Sikorski, PLoS One, 2009, 4, e5053.
33. S. R. Shaikh, V. Cherezov, M. Caffrey, S. P. Soni, D. LoCascio, W. Stillwell and S. R. Wassall, J. Am. Chem. Soc., 2006, 128, 5375–5383.
34. S. J. Marrink, A. H. de Vries and A. E. Mark, J. Phys. Chem. B, 2003, 108, 750–760.
35. T. Carpenter, P. J. Bond, S. Khalid and M. S. P. Sansom, Biophys. J., 2008, 95, 3790–3801.
36. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
37. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
38. P. Dynarowicz-Łątka and K. Hąc-Wydro, Colloids Surf., B, 2004, 37, 21–25.
39. S.-S. Qin, Z.-W. Yu and Y.-X. Yu, J. Phys. Chem. B, 2009, 113, 8114–8123.
40. S. J. Marrink, J. Risselada and A. E. Mark, Chem. Phys. Lipids, 2005, 135, 223–244.
41. S. Sonnino, A. Prinetti, L. Mauri, V. Chigorno and G. Tettamanti, Chem. Rev., 2006, 106, 2111–2125.
42. W. J. Allen, J. A. Lemkul and D. R. Bevan, J. Comput. Chem., 2009, 30, 1952–1958.
43. J. Wohlert and O. Edholm, J. Chem. Phys., 2006, 125, 204703.
44. D. H. de Jong, C. A. Lopez and S. J. Marrink, Faraday Discuss., 2013, 161, 347–363.
45. H. Koldsø, D. Shorthouse, J. Hélie and M. S. P. Sansom, PLoS Comput. Biol., 2014, 10, e1003911.
46. P. Balgavý, M. Dubničková, N. Kučerka, M. A. Kiselev, S. P. Yaradaikin and D. Uhríková, Biochim. Biophys. Acta, Biomembr., 2001, 1512, 40–52.
47. S. Tristram-Nagle, R. Zhang, R. M. Suter, C. R. Worthington, W. J. Sun and J. F. Nagle, Biophys. J., 1993, 64, 1097–1109.
48. F. de Meyer and B. Smit, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 3654–3658.
49. H. A. Rinia, M. M. E. Snel, J. P. J. M. van der Eerden and B. de Kruijff, FEBS Lett., 2001, 501, 92–96.
50. L. Janosi, Z. Li, J. F. Hancock and A. A. Gorfe, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 8097–8102.
51. S. May, Soft Matter, 2009, 5, 3148–3156.
52. C. Rosetti and C. Pastorino, J. Phys. Chem. B, 2012, 116, 3525–3537.
53. A. Filippov, G. Orädd and G. Lindblom, Biophys. J., 2004, 86, 891–896.
54. G. Lindblom and G. Orädd, Biochim. Biophys. Acta, Biomembr., 2009, 1788, 234–244.
55. O. Bakht, P. Pathak and E. London, Biophys. J., 2007, 93, 4307–4318.
56. W. Guo, V. Kurze, T. Huber, N. H. Afdhal, K. Beyer and J. A. Hamilton, Biophys. J., 2002, 83, 1465–1478.
57. J. M. Holopainen, A. J. Metso, J.-P. Mattila, A. Jutila and P. K. J. Kinnunen, Biophys. J., 2004, 86, 1510–1520.
58. J. Aittoniemi, P. S. Niemelä, M. T. Hyvönen, M. Karttunen and I. Vattulainen, Biophys. J., 2007, 92, 1125–1137.
59. A. G. Ostermeyer, B. T. Beckrich, K. A. Ivarson, K. E. Grove and D. A. Brown, J. Biol. Chem., 1999, 4, 34459–34466.
60. X. Xu and E. Londin, Biochemistry, 2000, 39, 843–849.
61. R. Schroeder, E. London and D. Brown, Proc. Natl. Acad. Sci. U. S. A., 1994, 91, 12130–12134.
62. E. London, Biochim. Biophys. Acta, Mol. Cell Res., 2005, 1746, 203–220.

---

Footnote

† Electronic supplementary information (ESI) available: Fig. S1, individual area per lipid for Chol in the three systems; Fig. S2, membrane thickness; Fig. S3, logarithmic histogram of the Chol translocation times; Fig. S4, time series for the slowest Chol translocations; Fig. S5, time evolution of the diffusion coefficient for DSPC; Fig. S6, time evolution of the diffusion coefficient for POPE; Fig. S7, time evolution of the diffusion coefficient for Chol. See DOI: 10.1039/c5ra02196k

---