Long chain sphingomyelin depletes cholesterol from the cytoplasmic leaflet in asymmetric lipid membranes

The transbilayer distribution of cholesterol (CHL) in complex asymmetric lipid membranes remains controversial, with contrasting investigations suggesting that there is more CHL either in the exoplasmic, outer leaflet (OL) or the cytoplasmic, inner leaflet (IL) depending on cell type or model, membrane composition, and method of investigation. Here, we launch systematic coarse-grained molecular dynamics simulations to investigate the impact of the sphingomyelin (SM) acyl chain length upon CHL distribution in asymmetric lipid membrane mixtures which account for the variation of the most abundant headgroups and acyl chain unsaturation in the two membrane leaflets. We find that there is always more CHL in the OL, but longer chain SM depletes more CHL from the IL than short chain SM in simple membrane mixtures containing SM and 16 : 0, 18 : 1 phospholipids. The difference between longer and shorter chain SM is neutralised in a more complex asymmetric membrane, where there are more saturated tails in the outer leaflet. We propose that interdigitation of long-chain SM into the opposing IL pushes cytoplasmic CHL towards the OL, but higher chain saturation of the outer leaflet compensates for the effect of SM chain length.


Introduction
CHL is an important and the most abundant (30-50% of all lipids) molecule in the plasma membrane (PM) that surrounds all eukaryotic cells. It has several vital roles in the PM e.g. maintaining membrane integrity, modulating membrane uidity over the range of physiological temperatures, forming PM nanodomains which putatively serve as hotspots of biochemical reactions and binding events, 1 and inuencing the conformation and function of several key classes of membrane proteins including G-protein coupled receptors, mechanosensitive ion channels and ion pumps. [2][3][4][5][6] While the concentration of the CHL in the cell and different cellular membranes has been well studied, its trans-bilayer distribution remains controversial. 7 An even distribution of CHL molecules in the two leaets is expected in a at, symmetric bilayer. However, the PM of eukaryotic cells is highly asymmetric with respect to the lipid composition of the two leaets (trans-bilayer asymmetry). The asymmetry drives important cellular functions and pathways. Most of the sphingomyelin (SM) and phosphatidycholine (PC) lipids are within the outer, exoplasmic leaet (OL) while the inner, cytoplasmic leaet (IL) mostly consists of the phosphatidylethanolamines (PE), phosphatidylserine (PS) and phosphatidylinositol (PI) lipids. 7,8 Trans-bilayer asymmetry is not an equilibrium conguration of a membrane (at equilibrium, an even distribution is expected). The cell uses active processes carry out by ippase proteins to maintain trans-bilayer asymmetry because there is a large free energy barrier associated with the translocation of a charged lipid from one bilayer leaet to the other across the membrane hydrophobic core. The transbilayer distribution of the major phospholipid components of the PM can now be carefully determined using highly sensitive lipidomics experiments. 9 However, data pertaining to the distribution of CHL across the two bilayer leaets remains contradictory owing to the low free energy barrier of CHL to ip op across the two leaets. Some results show that most of the CHL molecules are concentrated in the IL while others suggest that the majority exist within the OL, reviewed in ref. 10 and 11. It has also been suggested that CHL distribution depends on the cell type and age. 12 The distribution of CHL across the bilayer leaets is likely to be determined by a passive thermodynamic equilibrium distribution that equalizes the CHL chemical potential in the two leaets. An active process that maintains CHL asymmetry (if it exists) is unlikely considering the rapid ip-op of the CHL, and remains to be discovered. Based on the above arguments, a theoretical model suggested that most CHL molecules accumulate in the IL due to the negative curvature of the PE lipids present in this leaet. 10,11 Another mechanism that can also drive CHL to the IL is the presence of the PS and PI lipids, because a coarse grained (CG) simulation study demonstrated that CHL has a higher affinity to the charged leaets. 13 However, the same simulation study suggested that CHL prefers saturated lipids, which are more abundant in the OL. Another factor that can drive CHL to the OL is the presence of SM lipids that exclusively exist in this leaet. Theoretical models predict that 80% of CHL resides in the OL owing to the presence of SM, but contributions from the bending energy reduce this abundance to 63%. 14 The majority of SM in mammalian cells in the outer leaet have a long C24 acyl chain, 9 which can interdigitate into the opposing leaet. 15,16 However, there are also more than 4% short-chain SM in the outer leaet. In the inner leaet, the small fraction (less than 0.1%) of SM are all exclusively short-chain lipids. 9 Mutations which replace C24 SM by C18 SM have physiological consequences related to obesity and insulin resistance. 17, 18 We previously showed that a delicate balance between SM-CHL hydrogen bonding and interactions with the solvent result in weaker interactions with C24 SM with CHL, compared to shorter chain C18 SM lipids. 16 At the same time, free energy calculations showed that CHL preferred the OL in the presence of C18 SM and the IL in the presence of C24 SM. 16 However, no PE or PS lipids were present in the lower leaet in these all-atom simulations. In addition, all-atom simulations suffer from the constraint that CHL cannot spontaneously ip op on the simulated timescales (100 s of nanoseconds). 19 Overall, one can argue that SM and saturated lipids pull more CHL into the OL, while PE and charged lipids attract more lipids to the IL. The competition between these determinants decides the equilibrium distribution of the CHL in the PM. Here, we exploit the possibility of spontaneous CHL ip-op in molecular dynamics (MD) simulations employing the CG MARTINI force eld for lipids, to investigate the effect of SM chain length on the distribution of CHL in the two bilayer leaets with compositions reminiscent of biological membranes, but simplied to eliminate minor lipid species. The ip-op rate of CHL in simulations of CG lipid bilayers can range from 1-20 ms À1 , 20-23 depending on the presence of unsaturated lipids (more unsaturated lipids increase ip-op rates), cholesterol content (more cholesterol lowers ip-op rates) and overall bilayer composition (overall, a more ordered bilayer will reduce CHL diffusion and ip-op rates). Simulations of asymmetric membranes have to be set up with care to ensure either that the tension in each leaet is zero, or equivalently, the area of each leaet is identical. 24 When investigating CHL asymmetry, if the overall area of each leaet is not calculated carefully, CHL will simply accumulate into the leaet containing too few lipids. We therefore employ a careful methodology which minimizes this area difference (described in detail in the Methods section). We nd that long-chain SM lipids interdigitate into the IL, and push more CHL into the OL to achieve mechanical equilibrium and balance the free volume in the two leaets.

Composition of asymmetric bilayers
We chose our bilayer compositions based on the lipidome of human erythrocytes, 9 simplied to eliminate minor lipid species and all proteins. Two compositions were chosen. First, a simple asymmetric membrane was modelled with predominantly PS and PE in the IL, and PC and SM predominantly in the OL. The lipid acyl chains were set to 16 : 0-18 : 1 for PS, PE and PC, while the SM was either all 18 : 1-18 : 0 (DPSM) or 18 : 1-24 : 0 (DSSM). The overall CHL fraction was 0.29, similar to prior simulations of plasma membrane models. 25 The nal system compositions are shown in Table 1. To account for complexities in acyl chain length, a complex asymmetric bilayer was also setup with the composition shown in Table 1. Lipids present at less than 5% were not included.

Construction of asymmetric bilayers
Unlike symmetric bilayers, the construction of asymmetric MD models of bilayers is challenging because the total area of the two leaets must be as close as possible to each other to avoid artifacts arising from imperfect lipid packing. For example, if one leaet has too few lipids, cholesterol will accumulate in this leaet to ll up the free volume and establish mechanical equilibrium. To avoid such artifacts, we devised the following workow to eliminate area mismatch as much as possible.
To calculate optimal area per lipid (APL) of each component, we simulated symmetric bilayers, where each monolayer had the desired composition of a single leaet of the target asymmetric bilayer. For example, for the simple systems, we ran simulations of symmetric POPC-SM-CHL bilayers of composition of the OL to obtain the most likely area per lipids of POPC, SM and CHL, and a symmetric POPC-POPE-POPS-CHL bilayer of composition of the IL, to obtain the most likely APL of POPC, POPE and POPS in the IL. The simulation box size was set to 14 Â 14 Â 10 nm for these simulations, and APLs were obtained from Voronoi triangulation using apl VORO. 26 From these APLs, the number of lipids required to model a bilayer of a target total area was obtained. We then constructed an asymmetric bilayer using the number of lipids calculated in the previous step. Such an asymmetric bilayer should have minimal area mismatch between the two leaets, presuming that the CHL content is identical in the two leaets. The preceding assumption is certainly not correct if CHL distribution is asymmetric, and we account for this in the following manner.
Aer a 10 ms production run, a trans-bilayer CHL asymmetry is created by overall ip-op of CHL molecules to equalize the chemical potential of CHL in the two leaets. However, such an asymmetry agains creates an area mismatch of DA chol in the two leaets, because CHL occupies additional area in the leaet with more CHL.
To eliminate area mismatch again, we ll the CHL-depleted leaet with additional lipids which occupy an area corresponding to DA chol , and repeat another 10 ms production iteration. The overall CHL concentration is kept constant across all iterations. Aer several such iterations, it is expected that the CHL distribution will stabilise across the two leaets.
A second approach to minimising the area mismatch is to choose different initial ratio of CHL in the two leaets in the starting conguration of the simulations. We launched simulations where the ratio between the number of CHL in the OL to the total number of CHL in the IL varied from 0.5 to 0.62, i.e. CHL was always richer in the OL, as indicated by the rst set of simulations (see Results section). Note that the number of phospholipids for each leaet remained the same. For each setup, we simulated 5 replicas with different initial velocity distributions. Only the systems with simple lipid compositions were simulated in this way. This set of simulation is referred to as the REPLICA simulations hereaer.
Overall, we have ran close to 770 ms of simulations in MARTINI time, corresponding to a real sampling time of close to 3.1 ms. 27

Simulation details
The coarse-grained (CG) MARTINI force eld 27-29 is used to model the lipid bilayers (required les were obtained from the MARTINI web site). 30 In MARTINI, 4 heavy atoms are typically represented by a single particle. For example, CHL is modelled by 8 and POPC by 12 beads. A CG force eld is necessary to permit CHL ip op to measure trans-bilayer CHL asymmetry. All lipid bilayers are modeled using the insane.py script obtained from the MARTINI web site, 30 in combination with home-made Matlab scripts. The bilayers are solvated with nonpolarizable water beads along with antifreeze beads. 27 Na + and Cl À ions are present to ensure a salt concentration of 150 mM. All simulations are performed using GROMACS 2019.x. 31, 32 The system was equilibrated using the Berendsen coupling 33 for both temperature and pressure. Production runs are carried out in the NPT ensemble for 2 ms using a time step of 20 fs. The temperature is controlled using the Nosé-Hoover thermostat, 34,35 with a reference value of 310 K, and the pressure is coupled semi-isotropically using the Parrinello-Rahman barostat, 36 with a reference value of 1 bar. The van der Waal's forces are shied to zero from 0.9 nm to 1.2 nm. Electrostatic interactions are calculated using the reaction-eld method as recommended. 37 The CHL distribution is measured from the z-coordinate (bilayer normal) of the hydrophilic ROH bead of CHL. Note that CHL can occupy the bilayer center in MARTINI simulations, but for the sake of simplicity, all CHL molecules are assigned either to one leaet or the other based on the position of the ROH bead. No CHL molecules are classied as being in the center of the bilayer. Lipid tail order parameters are calculated from a modied version of scripts available on the MARTINI web pages. 30 All other quantities are measured using GROMACS tools. All reported properties are calculated for the last ms of simulation.

CHL distribution
For the simple systems, the presence of SM in the OL attracts a higher fraction of CHL in this leaet compared to the IL (Fig. 1). The distribution ratio of CHL in the OL versus the IL converges aer 9 and 12 iterations of 10 ms each for the simple DSSM and DPSM systems respectively. More CHL partitions into the OL in the presence of longer chain SM (C24), than when C18 SM is present instead (Fig. 1). Since CHL ip ops rapidly in MARTINI simulations, and because CHL molecules close to the bilayer center are assigned to either leaet, there are some uncertainties associated with the data shown in (Fig. 1). However, it is clear that rst, more CHL partitions into the OL, and secondly, that SM chain length has a signicant impact on the CHL distribution ratio. The ip-op rates in our simulations are 1.95 and 3.65 per ms for the DPSM and DSSM simple simulations respectively. These numbers are averaged over all CHL molecules and correspond to a large number of ip-op events. The higher rate in the DSSM simulation is most likely the result of more CHL ip-opping into the outer leaet. The  Table 1. The yaxes show the ratio between the number of CHL molecules in the OL to the total number of CHL molecules in the system. If the ratio is 0.5, then CHL is evenly distributed in the two leaflets. In both cases, there is more CHL in the OL, than in the IL. The error bars are smaller than the size of the symbols and are therefore invisible.
numbers compare favourably to previous estimates (5.48 per ms) for systems without SM lipids. 20,38 With SM, the ip-op rate is expected to be slower owing to further ordering of the membrane.
For the complex systems, the distribution ratio stabilises aer 9 iterations, but the trend is different from the simple systems (Fig. 1). For DSSM, the distribution ratio is similar for the simple and complex systems. For DPSM, however, more CHL distributes in the OL in the complex system. One factor which leads to more CHL in the OL is that there are less saturated lipid acyl chains in the IL ($42% of phospholipid tails) in the simple systems, compared to the complex systems ($75% of phospholipid tails), and CHL is known to prefer saturated tails over unsaturated tails. 13 Transbilayer asymmetry in the distribution of saturated tails led to a signicantly increased proportion of CHL in the leaet containing saturated chains in prior CG simulations. 21,39 The increased saturation of the OL lipids scrambles the difference between DPSM and DSSM in the complex systems.
With the exception of two points where the difference is not statistically relevant, the differences between DPSM and DSSM are also retained for the REPLICA systems where different initial CHL ratio were enforced in the starting simulated congurations (Fig. 2). In principle, different initial CHL ratio should lead to the same nal CHL ratio aer sufficiently long sampling. However, different initial percentages of CHL in the two leaets lead can cause different initial areas of the two leaets and lead to different nal CHL ratio. It is therefore paramount that the simulations are carried out carefully while minimising the initial area mismatch as much as possible, as is done in the simulations with several consecutive iterations (data in Fig. 1) when the area of the OL and IL are kept equal at each instance. Furthermore, since there is always more SM than CHL in the OL, it is possible that a higher number of initial CHL in the OL leads to the formation of more stoichiometric SM-CHL complexes, thus retaining more CHL in the OL towards the end of the simulation. It is worth reiterating that the apparent lack of convergence of the data in Fig. 2 is not a consequence of limited sampling, but of problems arising out of initial area mismatch in the two leaets.
Transbilayer asymmetry in acyl chain saturation is identical in all our simple systems. The presence of SM in the OL, and the length of the SM chain are the other two key driving forces which determine the difference in the CHL distribution ratio. In the remainder of this article, we will analyse the simple systems further in an attempt to converge on a hypothesis of why longchain SM drives more CHL into the OL, in apparent contrast to our earlier ndings using all-atom simulations in conjunction with a host of analytical techniques. 16 For this, we rst investigate if the structural properties of the bilayer are different in the DPSM (C18) and the DSSM (C24) simple membrane simulations.
For analysing bilayer properties, we chose the nal iterations of the DPSM and DSSM simple system simulations. For timeaveraged properties, the last 1 ms of the simulations were analysed.

Bilayer properties
The thickness of the bilayers is calculated from the average distances between the phosphate beads in the two leaets. The thickness of the longer chain DSSM containing bilayers is expectedly higher, by about 0.5 nm (Fig. 3). There is an increase in the thickness over the rst few iterations, probably resulting from more CHL moving the OL, and forming nanodomains with SM lipids (see Fig. 4 later), straightening them in the process, resulting in increased membrane thickness.
The lipid tail order parameters demonstrate expected trends (Fig. 3B). The SM tails have higher order parameters compared to the POPC lipid tails. Slightly higher order parameters in the DSSM containing bilayers near the bilayer center result from the interdigitation of the DSSM acyl tails into the opposing IL (see next section), resulting in more ordered bilayers near the bilayer center. Interdigitation and the resulting higher order parameters also contribute to the increased thickness of the DSSMcontaining bilayers. However, neither the thickness or the  CHL and SM form stable nanodomains in both the DPSM and DSSM systems, as apparent from 2D-density plots (Fig. 4C) and (Fig. 4D). Although phase separation is not complete, several nanodomains 2-3 nm in size are visible for both systems (Fig. 4). The pairing of the CHL headgroup and the SM headgroup is also apparent form radial distribution functions between the ROH bead of CHL and the sphingosine backbone bead AM1 of SM (Fig. 4A) and (Fig. 4B). The interactions of CHL with the AM1 bead of SM are comparable with the POPC glycerol backbone GL1 bead. The binding of CHL to PS lipids is weaker than the binding to either the POPC or SM head groups.
The C24 DSSM acyl chain protrudes into the opposing IL, while the C18 DPSM does not (Fig. 5). Since headgroup interactions, domain formation, thickness and order parameters are identical in the DPSM and DSSM systems. Furthermore, PS lipids, which are known to attract CHL to the IL 13 are present in equal amounts in both the DPSM and DSSM systems. Therefore, we conclude that the interdigitation of lipid tails in the IL makes the IL denser, and reduces the free volume available for CHL. CHL is the only molecule which can spontaneously ip-op in the simulations on the ms timescale, and therefore CHL accumulates in the OL when the C24 SM tail interdigitates into the IL. The hypothesis is in conict with our previous ndings, which were based on free energy calculations from all-atom simulations. 16 The lipid composition in the said systems did not contain PS or PE lipids and therefore dissimilar from the ones we have used in the current CG simulations. More CHL in the OL (54%) was also observed in CG simulations of a bilayer with composition mimicking a plasma membrane. 25 There are several key differences in the said simulations from our current work. Firstly, we have focussed on the differences between the effects of long-chain and short-chain SMs, while the systems in 25 more resemble a realistic plasma membrane containing a mixture of short-chain and long-chain SM. Secondly, our systems are less likely to suffer from an initial area mismatch owing to reasons discussed earlier. Finally, our systems are much simpler in composition, with all lipids with an abundance of less than 5% being ignored.
The protein content in a real PM is nearly 30%, and CHL affinity to peripheral of specic regions of TM proteins can inuence the overall distribution of CHL. Although inclusion of proteins could not be performed in the current set of simulations, our scheme of carefully equilibrating the systems to keep the areas of each monolayer constant can be used in the future to analyse complex membrane systems also containing proteins.

Conclusions
Overall, our data indicate that long-chain SM in the OL attracts more CHL to the OL, than short-chain SM does. There are, however, many other compositional determinants which will ultimately determine the CHL distribution ratio in both experiments and simulations. Since lipid composition has such an important role to play, the distribution of CHL will be sensitive to the cell and organelle type, as has been proposed in the past. Our simulations indicate that short-chain SM, when present in cells in diseased conditions, result in less CHL accumulation in the OL. A different CHL distribution can in principle affect nanodomain formation in complex membrane environments, and also inuence signaling processes associated with these domains.

Conflicts of interest
There are no conicts to declare.   5 The density of the terminal SM bead along the bilayer normal for the DPSM and DSSM systems. The arrow highlights the tail in the distribution in the DSSM system resulting from the interdigitation of this bead into the opposing IL.