Hydrophobic Compounds Reshape Membrane Domains

Cell membranes have a complex lateral organization featuring domains with distinct composition, also known as rafts, which play an essential role in cellular processes such as signal transduction and protein trafficking. In vivo, perturbations of membrane domains (e.g., by drugs or lipophilic compounds) have major effects on the activity of raft-associated proteins and on signaling pathways, but they are difficult to characterize because of the small size of the domains, typically below optical resolution. Model membranes, instead, can show macroscopic phase separation between liquid-ordered and liquid-disordered domains, and they are often used to investigate the driving forces of membrane lateral organization. Studies in model membranes have shown that some lipophilic compounds perturb membrane domains, but it is not clear which chemical and physical properties determine domain perturbation. The mechanisms of domain stabilization and destabilization are also unknown. Here we describe the effect of six simple hydrophobic compounds on the lateral organization of phase-separated model membranes consisting of saturated and unsaturated phospholipids and cholesterol. Using molecular simulations, we identify two groups of molecules with distinct behavior: aliphatic compounds promote lipid mixing by distributing at the interface between liquid-ordered and liquid-disordered domains; aromatic compounds, instead, stabilize phase separation by partitioning into liquid-disordered domains and excluding cholesterol from the disordered domains. We predict that relatively small concentrations of hydrophobic species can have a broad impact on domain stability in model systems, which suggests possible mechanisms of action for hydrophobic compounds in vivo.


Introduction
Biological membranes are both chemically and structurally heterogeneous. The constituent lipids can self-organize in domains [1], which differ in chemical composition and in physical properties, including structural, dynamic, and elastic properties. Domains have a functional role in cells: membrane proteins partition preferentially to one specific domain (or to domain boundaries) and carry out their function correctly only when in the appropriate environment -as expressed by the raft concept [1]. Membrane lateral organization is involved in biological processes such as membrane fusion [2,3], signal transduction [4], protein trafficking [5], and viral infection [6,7]. Alterations of the membrane lateral organization have been identified in pathologies like allergies and the Alzheimer disease [8], and have been linked to the mechanism of action of general anesthetics [9,10]. Understanding the determinants of domain stability in vivo is therefore of paramount importance in biomedical sciences. Yet, characterization of raft domains in vivo is challenging because of the small size of the domains, which are typically smaller than optical resolution [11]. In model systems (i.e., vesicles), instead, domains are usually larger and can even coalesce to yield macroscopic phase separation. For this reason, model systems are often used to study membrane lateral organization [11].
Among model systems, the most frequently used are ternary mixtures of cholesterol and two lipids with different melting temperatures, as they show liquid-ordered (L o ) -liquid-disordered (L d ) phase coexistence, similar to cell membranes [12].
Compounds with sufficiently high affinity for membranes can modulate biological function by virtue of membrane-mediated effects [13][14][15], including the alteration of membrane lateral organization. Recent studies have shown that, in model membranes, phase coexistence is affected by a variety of compounds. For instance, some lipids [16], vitamin E [17], and n-alcohols [10] destabilize phase separation in ternary lipid mixtures. On the contrary, transmembrane helical peptides [18], benzyl alcohol [17], and polystyrene [19] stabilize phase separation. It is unclear which chemical or physical properties of the solutes determine stabilization or destabilization of phase separation. Systematic studies on the effect of solutes on membrane lateral organization are lacking. Moreover, the mechanisms of stabilization and destabilization of domains are not understood.
In the present report, we describe the effect of different hydrophobic compounds on lipid mixing in phase-separated membranes. Hydrophobic compounds partition largely to the interior of lipid membranes, hence they do affect many membrane properties. Hydrophobic compounds are extremely common in commercial products and in the environment; for instance, they are used as fuels in combustion engines, as solvents in industrial processes, and as scaffolds in drugs. Also, they are building blocks for many industrial polymers and they are found in the atmosphere as pollutants (e.g., in products of incomplete combustion of fossil fuels). We determine the effect of hydrophobic compounds on phase separation using coarse-grained (CG) molecular dynamics (MD) simulations of L o -L d phase-separated lipid membranes. We focus on six different hydrophobic solutes covering a wide range of sizes and a variety of chemical structures: cyclohexane, octane, hexadecane, benzene, C 60 fullerene and polystyrene. All solutes partition to the interior of the membrane but, remarkably, they show very different lateral distributions. We identify two distinct groups with different lateral distributions, and we show that they have opposite effects on lipid mixing. Finally, we determine the mechanism of action for both groups of molecules.

Lipid mixing: Effect of aromatic vs. aliphatic solutes
We used the MARTINI coarse-grained (CG) force field [20,21] to simulate model membranes consisting of dipalmitoyl-phosphatidylcholine (DPPC), dilinoleyl-phosphatidylcholine (DLiPC), and cholesterol, at 42:28:30 molar ratio. At a temperature of 295 K, in the absence of solutes, the membrane showed phase separation into a liquid-ordered (L o ) domain, comprising mostly DPPC and cholesterol, and a liquid-disordered (L d ) domain, comprising mostly DLiPC, as reported previously [22]. Due to periodic boundary conditions used in the simulations, the domains organized in stripes along one axis of the box (Fig. 1). These stripes persisted during the simulation, yet the interfaces were dynamic and lipid molecules exchanged between the L o and L d phase.
We then carried out simulations of the same membrane in the presence of six different hydrophobic solutes: octane, hexadecane, cyclohexane, benzene, fullerene, and polystyrene. These solutes are chemically diverse, as they include linear alkanes of different size, a cyclic alkane, a common aromatic hydrocarbon, a wellknown carbon nanoparticle, and a common industrial polymer. For each solute, we performed two simulations at low solute concentration (3.3% solute/lipid molar ratio) and two simulations at high solute concentration; due to the very different molecular weight of the solutes, for the high concentration we chose to use a common solute/lipid mass ratio of 4.8%. For polystyrene, we only considered simulations at high concentration. Additional simulations where performed for some solutes (see Table 1 for the complete list of simulations). Partition of the six solutes to the interior of lipid membranes is thermodynamically highly favorable, as shown in several previous studies [19,[23][24][25], but the time scales for permeation depend largely on the size of the particles: small molecules penetrate within a few nanoseconds [24] while large polymer particles require microseconds or more [19]. Since our interest was not in the kinetics of permeation but in the effect of the solutes on the membrane, we decided to start all simulations with the solutes placed inside the membrane, homogeneously distributed in the membrane plane.
To quantify phase separation, we calculated the DLiPC-DPPC contact fraction, f mix , defined as the fraction of DLiPC-DPPC contacts over the total number of contacts of DLiPC with all phospholipids (therefore not including cholesterol; see Methods). The DLiPC-DPPC contact fraction will tend to 0 at complete phase separation and will reach 0.61 at ideal mixing (equaling the DPPC molar fraction with respect to phospholipids only). In the absence of solutes, f mix was 0.1360.004, indicating strong phase separation. Addition of a small concentration of hydrophobic solutes had a minor effect on the DLiPC-DPPC contact fraction (0.11,f mix ,0.15, depending on solute type; see Table 1). Yet two trends were distinguishable: octane, hexadecane, and cyclohexane caused an increase in lipid mixing, while benzene, fullerene, and polystyrene caused a slight decrease in lipid mixing. These trends were more evident at high solute concentration: f mix reached 0.25-0.32 with the first group of compounds, and decreased to 0.10-0.11 with the second group (see Table 1). Visual inspection of the trajectories showed significant mixing (although not ideal mixing) in the presence of octane, hexadecane or cyclohexane, while domains were clearly separated in the presence of benzene, fullerene, or polystyrene ( Fig. 1).
The demixing effect induced by benzene, fullerene, and polystyrene appeared weaker than the striking mixing effect of octane, hexadecane, and cyclohexane. This is because the reference membrane was already phase-separated at 295 K. To assess domain stabilization by benzene, fullerene, and polystyrene, we carried out simulations at higher temperature. At the temperature of 325 K the system without any solute was no longer phase-separated, with f mix = 0.3460.02 (Fig. 2). Remarkably, the membrane remained clearly phase-separated in the presence of benzene, fullerene, and polystyrene, and the increase in lipid mixing at higher temperature was minor (Fig. 2). The DLiPC-DPPC contact fraction was only 0.2060.004 in the presence of fullerene at 325 K, and even less with benzene and polystyrene (Table 1). In contrast, in the presence of cyclohexane, lipids were already rather mixed at 295 K (at high concentration) and they mixed more at 325 K.
In summary, we observe a strong effect of all hydrophobic molecules on the stability of domains in phase-separated membranes, and we identify two groups of compounds with opposite effects on domain stability. The only obvious chemical property common within each group appears to be aromaticity (or the lack of it): all aromatic compounds promote lipid demixing, while all aliphatic compounds promote lipid mixing. What are, then, the mechanisms leading to such different effects on phase separation?

Author Summary
Cell membranes consist of a variety of lipids and proteins with inhomogeneous lateral distribution, forming domains with distinct composition and properties. These domains play a fundamental role in a number of biological processes, and perturbing them can have important effects on cellular functions. Some chemicals with high affinity for lipid membranes perturb membrane domains, but the link between properties of the chemicals and domain perturbation is not understood. The mechanisms of domain perturbation are also not understood. In the present work we use molecular simulations of model membranes to understand the driving forces and the mechanisms of domain perturbation by different chemicals. We explore the effect of six hydrophobic compounds, all of them rather simple and common but with different size, shape, and properties. We find that all hydrophobic compounds alter the stability of domains, but not all of them in the same way. We identify two groups of compounds with opposite effects: aromatic compounds stabilize domains, while aliphatic compounds destabilize them. Simulations also allow us to visualize, for the first time, the mechanism of domain perturbation -which is very difficult to assess experimentally. Our findings on model membranes suggest possible mechanisms of action for hydrophobic chemicals in living cells.

Aliphatic solutes favor lipid mixing by acting as linactants
To understand the driving forces for solute-mediated alterations of membrane lateral organization, we analyzed the spatial distribution of each solute within the membrane by calculating the solute-DLiPC contact fraction (Fig. 3), defined as the number of contacts the solute makes with DLiPC over the number of contacts it makes with all phospholipids (see Methods). If the solute is ideally mixed, the solute-DLiPC contact fraction will be equal to the molar fraction of DLiPC with respect to all phospholipids, i.e., 0.39; lower values indicate that the solute makes contacts preferentially with DPPC, while higher values indicate that the solute makes contacts preferentially with DLiPC. As shown above, high concentrations of aliphatic solutes caused lipid mixing; the solute-DLiPC contact fraction in those systems was close to 0.39, as expected -lipid mixing is concurrent with solute mixing. More interesting is the behavior at low solute concentration (at 295 K), when the domains remain phase-separated. At low concentration, aliphatic solutes showed a preference for DLiPC. However, this preference was small, with solute-DLiPC contact fractions close to 0.5, indicating that the solute made the same number of contacts with DPPC and DLiPC. Such situation occurs if the solute is either found in both phases (with a mild preference for DLiPC), or if it lies at the interface between them. Analysis of the density landscapes of the different components indicated that aliphatic solutes distribute preferentially at the interface between the domains. The preference was very clear for hexadecane, and it was observable also for octane and cyclohexane (Fig. 3). Such preferential distribution suggests that aliphatic compounds act as linactants.
To understand the mechanism of action of aliphatic compounds in more detail, we plotted the DLiPC-DPPC contact fraction vs. solute-DLiPC contact fraction as a function of simulation time. Figure 4A shows the results for one simulation with octane at high concentration at 295 K. Initially, by design, octane was evenly Once octane molecules reached the interface, they remained there while the lipids started mixing (2R3). After about 300 ns, the domains started to become blurry (3). While phase separation disappeared, octane mixed as well (3R4). Finally, both the lipids and the solute were mostly mixed (4). A very similar behavior was observed also with hexadecane and cyclohexane, although the detailed kinetics was different (Fig. S1). The sequence of events indicates clearly that all aliphatic compounds act as linactants, first moving towards the L d -L o interface and then destabilizing phase separation.
While the general mechanism of action was similar for all aliphatic compounds (Fig. S1), the kinetics of mixing depended on the nature and on the concentration of the solute: cyclohexane induced mixing faster than octane and hexadecane; also, in the presence of cyclohexane the transition (2R3) started before the solute-DLiPC contact fraction reached 0.5, i.e., before all the solute reached the interface. On the contrary, the time scale for mixing was longer in simulations with hexadecane; compared to octane and cyclohexane, hexadecane showed a higher affinity for unsaturated lipids.
Despite differences in the kinetics of lipid mixing, the extent of lipid mixing was remarkably similar for all aliphatic solutes at all concentrations, once concentrations were expressed as molar fractions (Fig. 4B). This indicates that the chemical potential of each lipid in the L o and L d phase did not depend on the type of aliphatic solute. In other words, the thermodynamics of lipid mixing was surprisingly independent of the nature of the solute.

Aromatic compounds favor phase separation by redistributing cholesterol
In contrast to aliphatic compounds, aromatic solutes such as benzene, fullerene, and polystyrene, stabilized the L o -L d phase separation. How did aromatics stabilize phase separation? To understand the underlying mechanism, we analyzed solute distribution in the membrane. Solute-DLiPC contact fractions for all aromatic compounds were close to 1, indicating a strong preference for the L d phase, as also confirmed by density landscapes (Fig. 5). An obvious potential mechanism to promote phase separation involves changes in the properties of the L d phase, where aromatics lie. For example, thinning of the L d phase would lead to an increase in thickness mismatch between the L d and L o domains, favoring phase separation. However, we found that all aromatic solutes caused an increase in the thickness of the L d domain (Fig. S3). As a result, the thickness mismatch between the two phases was actually reduced by these solutes, not increased. The largest reduction in thickness mismatch was observed in case of polystyrene, amounting to about 0.2 nm. Clearly changes in the thickness of the L d phase cannot explain the effect of aromatic compounds.
An alternative hypothesis is that aromatic solutes compete with the (small) fraction of cholesterol that resides in the L d phase. Visual inspection of the trajectories suggested that aromatic solutes replaced cholesterol in the L d phase (Fig. 5). Cholesterol-DLiPC contact fraction showed that few cholesterol molecules partitioned to the L d phase, both with and without added solutes. Yet, in the presence of aromatic compounds, the presence of cholesterol in the DLiPC-rich phase was significantly reduced, particularly at high temperature (see Table 1). We conclude that aromatic solutes, by partitioning into the L d domain, provide an additional driving force for cholesterol to enter the L o phase. As a result, the difference in order between the domains increases even further, and domain segregation becomes stronger.
Since the mechanism of action of aromatic compounds appeared to involve the displacement of cholesterol from the L d phase, we verified that this result does not depend strongly on the particular choice of the cholesterol-aromatic interaction. We carried out additional simulations with a modified force field, in which the strength of cholesterol-aromatic interaction was increased (see Methods for the details). We found that phase separation was about the same as with the original force field (Fig.  S4): in the presence of benzene, DLiPC-DPPC contact fraction was 0.1760.01, very similar to the contact fraction obtained with the regular force field in the same conditions (0.16). Moreover, solute-DLiPC and cholesterol-DLiPC contact fractions calculated with the original and with the modified force field were very similar (Table 1). Overall, our results indicate that phase separation and the mechanism of its stabilization by aromatics are robust with respect to reasonable variations in cholesterolaromatic interaction.

Discussion
Membrane lateral organization has paramount importance in cellular processes such as signaling, protein trafficking, and viral infection. Perturbations of membrane lateral organization can affect a large number of processes vital to the cell. Changes in domain structure can be brought about by modifications in membrane composition or by the addition of molecules that dissolve in the membrane. The effect of a few small molecules (alcohols [9,10], surfactants [17], anesthetics [9,10]) on membrane lateral organization has been studied experimentally in model systems. It has been observed that some molecules stabilize domains, while others destabilize them, but results are sparse, so it has not been possible to pinpoint the chemical or physical properties determining stabilization and destabilization of domains. Moreover, little is known on the mechanisms of lipid domain reshaping. Both the thermodynamics and the mechanisms of domain reshaping are difficult to study in living cells because of the small size of the domains and their highly dynamic nature.
Here we studied the effect of a set of common hydrophobic molecules on the lateral organization of model lipid membranes, consisting of saturated and unsaturated phospholipids, and cholesterol. Such membranes display clear phase separation between L d and L o phases at room temperature, and lipid mixing at higher temperatures -both experimentally and in MARTINI CG simulations. Our simulations predict that common hydrophobic compounds have major effects on lipid mixing in model membranes. Based on their effect on phase separation, the hydrophobic molecules selected for our study can be divided in two groups: (1) octane, hexadecane, and cyclohexane distribute preferentially at domain boundaries and destabilize phase separation; (2) benzene, fullerene, and polystyrene, instead, partition largely to the L d phase and stabilize phase separation. These predictions can be tested directly with experiments on model systems. Considering the diversity of the chemical structures used in our study, our conclusions are likely to be valid in a general way for purely hydrophobic compounds.  The persistence of phase separation in the presence of aromatic compounds at high temperature raises questions on the possibility that the systems might be trapped in metastable states. Based on the analysis of contact fractions, convergence requires about 1 ms in all simulated systems. Since our sampling is generally at least one order of magnitude longer, we expect phase separation to be well converged in all simulations. Yet, to guarantee that our simulations overcome potential metastable states, we repeated one simulation with benzene (high concentration, high temperature) starting from a well-mixed membrane (see Methods for details). Again, we observed phase separation within about 1 ms, and the contact fractions converged rapidly to the same values calculated in the original set of simulations (see Table 1). We conclude that persistence of phase separation in simulations with aromatic compounds is not due to limited sampling or the presence of metastable states.
Coarse-grained simulations provide both equilibrium and timedependent distributions of all species in a membrane; therefore they can be used to shed light on the mechanism of action of the different compounds -which is more difficult to access experimentally. For the first group of molecules (octane, hexadecane, and cyclohexane), we found that the mechanism of action is the one typical for linactants: those compounds tend to accumulate in the L d -L o interface region, which leads to a destabilization of the phase boundary [26]. For the second group (benzene, fullerene, and polystyrene), crowding of the L d phase prevents cholesterol from entering it, causing enrichment in cholesterol in the L o phase, particularly at high temperature. Cholesterol distribution in the L o phase has been associated to phase stabilization [27].
One of the goals of our study was to understand which chemical and physical properties of hydrophobic molecules determine their effect on domain stability. The compounds used in our study differ in several ways. Octane and hexadecane differ only in size, and they differ from the other compounds for the absence of ring structures. Hexadecane and cyclohexane are smaller than fullerene and polystyrene but bigger than benzene. Clearly the difference in domain remodeling behavior does not depend on the size of the solute. Nor it depends on cyclic nature of the compounds: both benzene and cyclohexane are cyclic (and of very similar size), but they have opposite effects on domain stabilization. Instead, the main discriminant between the two groups of compounds is aromaticity. The stronger affinity of aromatic compounds for the L d phase can be explained by p-p interactions between aromatic rings and double bonds in unsaturated acyl chains -which are captured by the force field in an effective way, through more attractive Lennard-Jones interactions. Aliphatic compounds, on the contrary, have higher affinity for saturated acyl chains (as expected based on experimental partitioning data [28]) but the dense packing of the L o phase prevents mixing of these solutes with the L o phase, as shown before for transmembrane peptides [29]. Dense packing of the L o phase appears to be responsible for preferential partitioning of aliphatic compounds at domain boundaries.
Together with the current study, there is a growing body of evidence indicating that small molecules can have a pronounced effect on lipid phase behavior. Like octane, hexadecane, and cyclohexane, also amphipathic molecules such as palmitoyloleoylphosphatidylcholine (POPC) [16] and vitamin E [17] distribute at the L o -L d interface and destabilize domains. Stabilization of phase separation has also been reported, for instance, in our previous work on polystyrene fragments of varying sizes [19], but also for transmembrane peptides [18] and for less hydrophobic solutes such as benzyl alcohol [17]. Exclusion of cholesterol from the L d phase due to crowding is likely to be the underlying mechanism of domain stabilization in all of these cases.

Conclusions
In conclusion, we showed that relatively small concentrations of six different hydrophobic compounds have a major impact on lipid domain stability in model membranes. Aliphatic compounds behave like linactants, accumulating at the interface between liquid-ordered and liquid-disordered domains and promoting lipid mixing, while aromatic compounds partition preferentially to liquid-disordered domains and stabilize phase separation. Both stabilization and destabilization of lipid domains can have an important impact on biological function. For example, it has been shown that, in vitro, raft-disrupting drugs can inhibit various cellular signaling pathways, including apoptitic pathways [30]. More studies are needed to understand how the complex interplay between lipids, proteins, and drugs affects signaling pathways in vivo. Nevertheless, our results on model systems shed light on the driving forces and the mechanisms of domain perturbation, and can be used to guide the rational design of drugs modulating phase separation. Knowledge of how hydrophobic molecules affect phase separation can also help understanding the side effects of drugs, and suggest possible mechanisms behind the toxicity of hydrophobic pollutants, such as hydrocarbons, air-borne carbon nanoparticles and nanoplastics.

System setup
We carried out all MD simulations at the coarse-grained (CG) level using the MARTINI force field [20,21,31]. The MARTINI force field is widely used for a large variety of membrane processes, including domain formation, as reviewed in refs [32] and [33]. We carried out simulations of lipid mixtures in water, containing 540 DLiPC, 828 DPPC, and 576 cholesterol molecules, as well as 21,880 water particles. The membrane was originally formed through self-assembly, by Risselada et al. [22] In the absence of solutes, at 295 K, the membranes yield phase separation and display a liquid-ordered (L o ) and a liquid-disordered (L d ) phase, which form stripes across the periodic box (Fig. 1).
We then simulated the same model membrane in the presence of six different hydrophobic solutes: octane, hexadecane, cyclohexane, benzene, C 60 fullerene, and polystyrene. We carried out the simulations at two solute concentrations. For the lower concentration, we used a constant solute:membrane molar ratio of 3.3% (i.e. 64 solute molecules). For the higher concentration, we used a constant mass ratio of 4.8% (based on real molecular masses). Except for the simulations with polystyrene, simulations started with the solute evenly distributed, on a grid, at the center of the membrane. The simulation time was between 6 and 30 ms (see Table 1). For polystyrene, we used simulations from a previous work by Rossi et al. [19]; in this case, the system contained 6 chains of 100 styrene residues each (PS100). PS100 chains formed compact clusters in water but dissolved once in the membrane interior, on a time scale of about 10 ms. In addition to the simulations above, some systems were simulated also at higher temperature or with additional solute concentrations. Simulations at higher temperatures were usually started from the same starting configurations used in simulations at lower temperature. One additional simulation was carried out in the presence of benzene starting with lipids completely mixed; in this case, the starting configuration was taken from the simulation with cyclohexane at 325 K, which showed a very high degree of mixing. Considering all simulations, the total sampling was over 680 ms. A list of all simulations performed is reported in Table 1.

Simulation parameters
The MARTINI [20,21] force field was used in all simulations. For simulations with fullerene, we used the fullerene model developed by Monticelli [24,34]. For simulations with polystyrene, we used the model by Rossi et al. [35]. One simulation was carried out with a modified force field, in which the strength of cholesterol-aromatic interaction was increased; namely, the SC4-SC1 interaction was increased from 3.1 kJ/mol to e = 3.5 kJ/ mol, while leaving all other interactions unchanged. Non-bonded interactions were calculated with a cut-off of 1.2 nm, on which we applied a shift function, starting at 0.9 nm for Van der Waals interactions and at 0 nm for Coulomb interactions. Charges were screened with a relative dielectric constant e rel = 15. A neighbor list (with a cut-off of 1.3 nm) was updated every 10 steps.
Simulations were run in the NPT ensemble. Pressure was coupled to 1 bar using a semi-isotropic barostat and the Parrinello-Rahman algorithm [36] (time constant of 4 ps and compressibility of 4.5610 25 bar 21 ). The temperature was coupled using the Bussi-Donadio-Parrinello thermostat [37] (time constant of 2 ps). We carried out most simulations at 295 K, and some additional ones at higher temperatures: 305 K, 315 K, and 325 K (see Table 1). We used the leapfrog integrator and an integration time step of 20 fs. The time step was reduced to 15 fs in simulations at temperatures of 315 K or higher, and 18 fs in all simulations with polystyrene. All simulations were carried out using the GRO-MACS software package (v4.5) [38].

Simulation analysis
Contact fraction. As a metric for phase separation, we used DLiPC-DPPC contact fractions, defined as: where c is the number of contacts between the two lipid species in subscript. Contacts were calculated only between PO4 beads of the lipids. We used a distance threshold of 1.1 nm, like in previous work by Domanski [18].
Solute and cholesterol lateral distribution were quantified by calculating solute-DLiPC and cholesterol-DLiPC contact fraction, respectively. These contact fractions are defined as: where X is either solute or cholesterol, and c is the number of contacts between the species in subscript. Contacts between cholesterol and DLiPC were calculated using only the PO4 bead of lipids and the ROH bead of cholesterol, with a distance threshold of 1.1 nm. For the solute-DLiPC contact fraction, we used all particles of the lipids and the solute, and a distance threshold of 0.8 nm. Averaging was done over the last 5 ms of each simulation, and errors were estimated by block averaging as implemented in GROMACS [38].
Density landscapes. We used two kinds of density landscapes to visualize the density of the different molecules in the plane of the membrane: the partial density landscape, and the DLiPC density fraction landscape. The partial density landscape was defined as the density of a given molecule calculated on a grid placed in the plane of the membrane (XY plane). The X and Y dimensions were divided in 50 bins each so the grid cells were about 0.460.4 nm. We averaged densities over the last 0.5 ms of the simulations. The DLiPC density fraction was defined as the fraction of DLiPC density over the total density of PC lipids for each cell; therefore, it can assume values between 0 (DPPC is the only lipid in that cell) and 1 (DLiPC is the only lipid type in that cell). The main interfaces are located where the DLiPC density fraction is 0.5. Landscapes were calculated using an in-house software freely available from our website (http://perso.ibcp.fr/ luca.monticelli, see also ref [39]).
Thickness calculations. Membrane thickness was calculated as the distance between the average positions of PO4 beads of the two leaflets on a grid. The X and Y dimensions were divided in 50 bins each, so the grid cells were about 0.460.4 nm 2 . For each cell, the thickness was averaged over the last 0.5 ms for each trajectory. The thickness of a phase in simulations with stripe domains was defined as the most frequent local thickness in the thickness landscape. Thicknesses were calculated with an in-house software freely available on our website (http://perso.ibcp.fr/luca. monticelli/). Left panels: histograms of membrane thickness in the absence and in the presence of different solutes. Thicknesses are calculated separately for the DLiPC and the DPPC components, based on the distance (along the bilayer normal) between phosphate groups in each leaflet. Line colors are the same as in the right panels. Right panels: difference in thickness between the DLiPC-rich and the DPPC-rich phases, in the absence and in the presence of different solutes. Numbers in parentheses indicate different replicas of the simulations. (TIFF) Figure S4 Robustness of the results. Snapshots of a single leaflet from simulations of systems with high concentration of benzene at 325 K, carried out with the original MARTINI force field (left panel) and the modified force field (right panel). Colors are the same as in Fig. 1: DPPC is colored in blue and DLiPC in red.