Molecular simulation of rapid translocation of cholesterol, diacylglycerol, and ceramide in model raft and nonraft membranes.

The translocation of lipids across membranes (flip-flop) is an important biological process. Slow exchange on a physiological timescale allows the creation of asymmetric distributions of lipids across cellular membranes. The location of lipids and their rate of exchange have important biological consequences, especially for lipids involved in cellular signaling. We investigated the translocation of cholesterol, ceramide, and diacylglycerol in two model bilayers using molecular dynamics simulations. We estimate half times for flip-flop for cholesterol, diacylglycerol, and ceramide of 20 μs, 30 μs, and 10 ms in a POPC bilayer, compared with approximately 30 min, 30 ms, and 30 s in a model raft bilayer (1:1:1 PSM, POPC, and cholesterol). Cholesterol has a large (54 kJ/mol) free energy of exchange between the POPC and raft bilayer, and therefore, it strongly prefers the more ordered and rigid raft bilayer over the more liquid POPC bilayer. Ceramide and diacylglycerol have relatively small free energies of exchange, suggesting nearly equal preference for both bilayers. This unexpected result may have implications for ceramide and diacylglycerol signaling and membrane localization.

of magnitude faster fl ip-fl op for cholesterol, ceramide, and diacylglycerol in the raft bilayer (hereafter referred to as "Raft") compared with a POPC bilayer. Cholesterol has a large free energy of exchange between the POPC and Raft bilayers, whereas ceramide and diacylglycerol have modest free energies of exchange or similar affi nities for both bilayers.

EXPERIMENTAL PROCEDURES
We simulated two model bilayers: pure POPC (64 POPC) and a 1:1:1 mixture of PSM, cholesterol, and POPC (22 lipids each) as a model Raft bilayer. Berger parameters were used for POPC ( 36 ). Parameters for PSM ( 37 ) and cholesterol ( 38 ) have been previously described. Palmitoyl-ceramide (CER) and palmitoyloleoyl-diacylglycerol (DAG) were built from PSM and POPC by replacing the headgroup with a hydroxyl group, with parameters from the GROMOS87 force fi eld ( 39 ). For water, we used the SPC model ( 40 ). Simulations were run with GROMACS ( 41 ). Water bonds and angles were constrained using SETTLE ( 42 ), and all other bonds using LINCS ( 43 ). The time step was 2 fs with neighbor list updates every 10 steps. Temperature was maintained at 323 K using the V-rescale algorithm ( 44 ). Lateral and normal to the bilayer plane, pressure was kept at 1 bar, using the Berendsen barostat ( 45 ). Lennard-Jones and Coulomb real-space interactions were cut off at 0.9 nm, with smooth particle mesh Ewald ( 46,47 ) for long-range electrostatics, with a fourth-order spline and 0.12 nm grid spacing.
We used umbrella sampling ( 48 ) to determine PMFs for moving lipids from bulk water to the center of the two model bilayers. We restrained the hydroxyl of each lipid at positions spaced by 0.1 nm along the bilayer normal with a force constant of 3,000 kJ mol Ϫ 1 · nm Ϫ 2 , resulting in 41 or 51 parallel simulations (to 4 or 5 nm from the bilayer center). Starting from an equilibrium bilayer, we pulled the lipid to its respective position with a force constant of 500 kJ mol Ϫ 1 nm Ϫ 2 for 1 ns, followed by 5 ns of equilibration with a 3000 kJ mol Ϫ 1 nm Ϫ 2 force constant and 15-60 ns of production simulation per window. The PMF was constructed using the weighted histogram analysis method ( 49 ). To increase computational effi ciency, we ran two umbrella samples in each simulation, one in either leafl et, with the lipid spaced by at least 4 nm. We have shown previously that the PMF for a DPPC lipid is not signifi cantly affected by using one or two umbrellas in the same system ( 50 ). Fig. 1 shows the kinetic model we used to estimate the rate of fl ip-fl op across a lipid bilayer, with k f the rate to move from equilibrium to the center of the bilayer and k d the rate to move from the bilayer center back to equilibrium. In a previous work ( 32 ), we showed that this method for estimating the rate for cholesterol fl ip-fl op matched rates directly observed from coarsegrained simulations.
To obtain k d , we fi rst pulled a lipid into the center of a bilayer and harmonically restrained it there for 1 ns with a 3,000 kJ mol ble domains, or signaling platforms, with the creation of ceramide by sphingomyelinase, is an example of a nonspecifi c mechanism of signaling ( 13 ). The localization and dynamics of signaling lipids is crucial to their function, particularly with respect to the distribution across the bilayer leafl ets, as lipids on one leafl et are effectively isolated from the other side. Ceramide created by neutral sphingomyelinase on the cytoplasmic leafl et is involved with a different signaling pathway than is ceramide generated on the outer leafl et by acid sphingomyelinase ( 13 ).
The translocation of lipids from one bilayer leafl et to another, or fl ip-fl op, is an important process for cellular growth, lipid transport, and signaling. It has been well established that the passive translocation of lipids with charged or zwitterionic headgroups is slow on a physiological timescale, with an estimated half time of hours to days (16)(17)(18). Proteins that facilitate lipid translocation, including fl ippases, fl oppases, and scramblases have been identifi ed, although identifi cation and functional characterization of membrane-perturbing integral membrane proteins is diffi cult ( 19 ). The much less polar headgroups of cholesterol, diacylglycerol, and ceramide suggest these lipids can translocate bilayers on a much faster timescale than phospholipids. Experimental evidence supports this hypothesis (20)(21)(22)(23)(24)(25)(26)(27), although the rates vary signifi cantly depending on the experimental conditions, often by orders of magnitude. These differences may be due to uncertainties with the interpretation of the complex experiments, but it seems more likely that they can be explained by a significant effect of differences in the structure and dynamics of the membranes themselves on the kinetics and thermodynamics of lipid fl ip-fl op.
Computer simulations of cholesterol-containing lipid bilayers are prevalent (28)(29)(30)(31). Previously, we investigated cholesterol fl ip-fl op and desorption using molecular dynamics (MD) simulations in pure phosphatidylcholine (PC) lipid bilayers ( 32 ). We found that cholesterol fl ips rapidly across PC bilayers, with a large gradient in rates between a polyunsaturated diarachidonylphosphatidylcholine (DAPC) bilayer, a pure DPPC bilayer, and a DPPC bilayer with 40 mol% cholesterol. A different study using the CHARMM force fi eld and 2D potentials of mean force (PMF) found similar free energy barriers to ours for cholesterol in DAPC, 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC), and DPPC bilayers ( 33 ). We showed previously that increasing cholesterol content in DPPC bilayers (0-40 mol%) dramatically increased the free energy barrier for DPPC fl ip-fl op and pore formation ( 34 ). Using similar methods, we also showed that saturated PC lipids with short tails (12 carbons) fl ip-fl op orders of magnitude faster than longer lipids (DPPC), through a pore-mediated mechanism ( 35 ).
Here, we have used atomistic MD simulations to investigate the interaction of cholesterol, palmitoyl-ceramide, and palmitoyloleoyl-diacylglycerol with an ordered and rigid raft-like bilayer [1:1:1 ratio of POPC, palmitoyl-sphingomyelin (PSM), and cholesterol] and with a pure POPC bilayer. We determined free energy profi les for moving a single lipid from the water phase to the center of the two different bilayers. From the profi les, we estimated orders We multiple the rate by one half to account for lipids that reach the bilayer center, then return back to their original leaflet. From the rate of fl ip-fl op, we then estimate the half time for fl ip-fl op using:

Equilibrium properties
Snapshots of the POPC and Raft bilayers during equilibrium simulations are shown in Fig. 2 . As expected, at equilibrium the headgroup of cholesterol, ceramide, and diacylglycerol resides at the water-lipid headgroup interface. The partial density profi les are shown in Fig. 2C . Both bilayers have peaks at the water-headgroup interface and troughs in the center of the bilayer. The Raft bilayer has a more complex density profi le, with multiple headgroup peaks due to multiple components. There is a trough at the center of the POPC bilayer and a deeper trough at the center of the Raft bilayer, likely because the two ordered leafl ets' tails are not able to interdigitate as much as more fl uid bilayers. The POPC bilayer is thinner and more fl uid than the Raft bilayer. Supplementary Fig. II shows order parameters for the palmitoyl tail of diacylglycerol in each bilayer. Supplementary Fig. II also shows that the Raft bilayer is much more ordered than the POPC bilayer. Fig. 3 shows PMFs for moving cholesterol, DAG, and CER from water to the center of POPC and Raft bilayers. All the PMFs share a number of common trends: a free energy minima at the equilibrium position of the headgroup in the bilayer and steep free energy slopes from equilibrium to bulk water or the bilayer From the free energy difference between equilibrium and the bilayer center and the rate for the lipid to move from the bilayer center to equilibrium, we can estimate a rate for fl ip-fl op. Using the equation : where k d is the rate to move from the bilayer center to equilibrium, ⌬ G center is the free energy difference, and k f is the rate to move from equilibrium to the bilayer center. We then determine the rate of fl ip-fl op (k fl ip ) using the rate to reach the bilayer center and the rate to move from the bilayer center to equilibrium:    Fig. 4 shows the number of hydrogen bonds formed between the cholesterol, diacylglycerol, or ceramide and the rest of the bilayer system as the lipid moves from water to the center of the bilayer. As the hydroxyl headgroup is moved into the bilayer core, water molecules move out of bulk solution into the hydrophobic interior to hydrogen bond with the lipid. There is an energetic balance between desolvating the fl ipping lipid and deforming the bilayer interface. We have observed similar behavior in the partitioning of charged and polar amino acid side chains in a dioleoylphosphatidylcholine (DOPC) bilayer ( 51 ). Moving ceramide from 0.4 nm ( Fig. 5B ) to 0.2 nm ( Fig. 5A ) from the POPC bilayer center causes water to stop interacting with the headgroup, illustrated by the average number of hydrogen bonds dropping to near-zero ( Fig. 4 ). The PMFs for the Raft bilayer have much steeper slopes compared with the POPC bilayer. The Raft bilayer's denser and stiffer structure makes it more expensive energetically to allow the polar lipid into its interior, thus center. The free energy for desorption (equilibrium to water) is much greater than the free energy barrier for fl ip-fl op. The important free energy differences are listed in Table 1 . The shaded regions in Fig. 3 represent d 2kt , the distance across the PMF at 2kT of energy, as a measure of the thermal fl uctuations accessible to each lipid normal to the plane of the bilayer. All three lipids have considerably larger fl uctuations in the POPC bilayer than in the Raft bilayer. Diacylglycerol has signifi cantly larger fl uctuations than either cholesterol or ceramide in both bilayers.

Flip-fl op
Both model bilayers are symmetric, and therefore, the free energy profi les shown in Fig. 3    of the bilayer with its tail in either leafl et, with slow exchange between the two conformations.
As ceramide and diacylglycerol are long and fl exible molecules, tilt angles are not well defi ned. Instead, we determined the average distance from the hydroxyl headgroup to the last carbon of the palmitoyl tail ( Fig. 6 ). For the POPC bilayer, there is a smooth decrease in this distance as the ceramide or diacylglycerol moves toward the bilayer center, refl ecting a heterogeneous distribution. In the Raft bilayer, both lipids have a deep trough in the head-to-tail distance at approximately 1 nm from the bilayer center, which suggests a compact conformation at this position. This distance corresponds to the free energy barrier for fl ip-fl op ( Fig. 3 ). At the center of the bilayer, the head-to-tail distance is close to the value at the equilibrium position, indicating mostly extended conformations. Fig. 5C shows the tail of ceramide interdigitated with one of the Raft bilayer leafl ets. Table 1    preventing translocation. At the free energy barriers, no hydrogen bonds between water and the fl ipping lipid are found, although diacylglycerol and ceramide do form one intramolecular hydrogen bond ( Fig. 4 ). At the center of the Raft bilayer, the numbers of hydrogen bonds for ceramide and cholesterol are not zero, due to an occasional water molecule reaching the center of the bilayer to hydrogen bond with the headgroup. This is likely possible due to the drop in density near the center of the Raft bilayer and could contribute to the small free energy troughs at the center of the Raft bilayer for cholesterol and ceramide ( Fig. 3 ). Fig. 5C shows ceramide hydrogen bonding with a water molecule at the center of the Raft bilayer.
The lipids change conformation as they move from water to equilibrium and then to the center of the bilayer. For cholesterol, we determined the tilt angle formed between its long axis and the normal to the plane of the bilayer ( Fig. 6 ) and the angle distribution ( Fig. 7 ). At equilibrium, cholesterol is oriented with its hydroxyl in water and its long axis perpendicular to the plane of the membrane. There is a slight increase in average tilt angle for the Raft bilayer compared with the POPC ( Fig. 6 ). As cholesterol is moved into the hydrophobic tail region, the angle distributions are poorly represented by an average value. Fig. 7A shows the angle distribution for cholesterol near the bilayer center, with a large peak at 150-180° for cholesterol at 0.0 nm and 0.2 nm from the bilayer center, but at 0-30° for 0.1 nm and 0.3 nm. There is also a small peak near 90° at 0.1 nm and 0.3 nm positions. Fig. 5D shows cholesterol at 0.2 nm from the center of the Raft bilayer, and facing the opposite leafl et at 0.4 nm from the bilayer center ( Fig. 5E ). Due to the order and rigidity of the Raft bilayer, cholesterol prefers to be perpendicular to the plane in the two model membranes. Using computer simulations of small bilayer patches allowed us to control the local concentration to simulate lipid fl ip-fl op in a Raft bilayer and a pure POPC bilayer . The Raft bilayer was both thicker and more ordered than the POPC bilayer. For all three lipids, the free energy barrier for fl ip-fl op increased considerably for the Raft bilayer compared with the non-Raft bilayer, and the rate for diffusion out of the interior decreased, resulting in orders of magnitude slower translocation across the thicker and less-fl uid Raft bilayer. All three lipids show orders of magnitude slower fl ip-fl op in the Raft bilayer compared with the non-Raft bilayer, with 8, 4, and 3 orders of magnitude difference for cholesterol, ceramide, and diacylglycerol. We fi nd that cholesterol has the slowest rate of fl ip-fl op in the Raft bilayer, compared with ceramide and diacylglycerol. In the pure POPC bilayer, the rate for fl ip-fl op follows the trend: cholesterol > diacylglycerol > ceramide.
Experimental values vary widely due to the use of different model systems, temperatures, labels, and techniques. For diacylglycerol and ceramide, studies also use different species with varied chain length and/or double bonds, which would affect the fl ip-fl op rate. Cholesterol fl ip-fl op is the most widely studied of the three, with experimental rate estimates ranging from hours ( 52 ) or minutes ( 27 ), to subseconds ( 26 ) or even submilliseconds ( 23 ). In a polyunsaturated DAPC bilayer, it was found that cholesterol is enriched in the bilayer interior, suggesting fast interleafl et diffusion ( 53 ).
Ceramide was found to fl ip rapidly across the plasma membrane of a live cell ( 54 ). The half time for unlabeled C-16 ceramide fl ip-fl op was estimated to be < 1 min in egg PC giant unilamellar vesicles ( 25 ). Subsequently, the same group showed that the half time for spin-labeled ceramide fl ip-fl op was < 1 min in both a liquid-disordered (egg PC) large unilamellar vesicle and a liquid-ordered (20:20:60 mol ratio of SM:POPC:CHOL) large unilamellar vesicle ( 24 ). This was in contrast to glucosylceramide and galactosylceramide, which both had signifi cantly reduced fl ip-fl op rates in bilayers with higher cholesterol concentrations ( 24 ). From a simple physiochemical perspective, it is hard to explain the discrepancy in fl ip-fl op behavior among these three types of ceramide. Our results suggest that ceramide fl ip-fl op in a liquid-disordered bilayer is faster than in a liquid-ordered bilayer, although we still fi nd a half time of 30 s in the liquid-ordered phase, which is faster than the experimental value ( 24 ).
In an erythrocyte membrane, 1,2-diacylglycerol was estimated to diffuse from the outer to inner leafl ets on the timescale of a minute ( 22 ). Fluorescently labeled diacylglycerol was found to equilibrate across DOPC vesicles on a fast timescale (minute), orders of magnitude faster than labeled phosphatidic acid ( 55 ). From NMR spectroscopy, the half time for 1,2-dilauroyl-glycerol was determined to be approximately 10 ms in egg PC vesicles ( 21 ). Overall, it is diffi cult to reconcile the varied experimental rates, and further work into studying lipid fl ip-fl op is warranted.
Our atomistic simulations illustrate the complexity for cholesterol, ceramide, and diacylglycerol translocation

Desorption
All three lipids have more hydrogen bonds in bulk water than at their equilibrium position in the bilayer ( Fig. 4 ). As the lipids move from equilibrium out into bulk water, their hydrophobic tails remain inserted in the bilayer to prevent exposure to water, resulting in a larger head-to-tail distance for ceramide and diacylglycerol. At a certain distance from the bilayer center, the tails stop interacting with the bilayer, and the lipid becomes fully desorbed. Ceramide and diacylglycerol show prominent drops in the head-to-tail distance as they stop interacting with the bilayer ( Fig. 6 ). For cholesterol, this is illustrated by the change in average angle from approximately 20° to 90° ( Fig. 6 ), the angle distributions for the Raft bilayer at 3.8 nm and 4.3 nm from the bilayer center ( Fig. 7B ), and the corresponding snapshots ( Fig. 5F, G ). The slope in the PMF for cholesterol desorption is much steeper for the Raft bilayer than the POPC bilayer, whereas ceramide and diacylglycerol have similar slopes for both bilayers. The free energy difference of a lipid at its equilibrium position and in bulk water is the free energy for desorption ( Table 1 ). Using the free energy in bulk water as reference point, we defi ne ⌬ G ex as the free energy for a single lipid to move from a Raft bilayer to a POPC bilayer; which is ⌬ G POPC -⌬ G Raft ( Fig. 1 ). Cholesterol has a strong preference for the Raft bilayer compared with the POPC bilayer, with a G ex of 54 kJ/mol. Diacylglycerol has a slight preference for the POPC bilayer compared with the Raft bilayer (G ex =-5 kJ/mol), whereas ceramide has a modest preference for the Raft bilayer (G ex = 7 kJ/mol).

Flip-fl op
We determined the free energy profi les for cholesterol, ceramide, and diacylglycerol translocation and desorption

Future directions
One of the sources of error in our PMF calculations are sampling problems, illustrated by the relatively large error bars for fl ip-fl op of the lipids in the Raft bilayer, in which lipid dynamics are slow. The simulation times are insufficient to sample the lateral diffusion of all lipids, but this potential problem is mitigated by the high concentrations of each of the lipid components and calculating two independent profi les for lipids or cholesterol in two different environments. There are also sampling problems as the lipid is moved out into bulk water ( 50 ). It is possible that the error bars underestimate the uncertainty in this region. The large jump in the head-to-tail distance ( Fig. 6 ) for ceramide and diacylglycerol between the tails inserted in the bilayer and the lipid folded up in bulk water, as well as the large error bars in this region illustrate the sampling problem. Supplementary Fig. III shows PMFs for cholesterol in the Raft bilayer using different lengths of simulation (from 5 ns to 60 ns). After around 20-30 ns, the PMFs converged to a reasonable statistical error for the purpose of this study, although orders of magnitude more sampling would allow much more accurate statistical analysis ( 58 ). Atomistic free energy profi les of lipids in bilayers are computationally demanding, and the present work is near the current state of the art in the fi eld. Increases in computational power, as well as improvements in free energy algorithms and enhanced sampling methods, will allow more accurate estimates to be performed on more complex systems.
Another potential source of error for our calculations, and simulations in general, are force fi eld inaccuracies, which are also diffi cult to assess. The parameters for ceramide, diacylglycerol, and cholesterol were derived in a consistent manner, with the same partial charges and atom types for the hydroxyl headgroup, which makes our comparison more reliable. The most important results we present are the differences in free energy, which likely mitigates potential force fi eld-related errors. We note that the free energy profi le we previously determined for cholesterol fl ip-fl op using the present parameters were similar to those calculated using the CHARMM force fi eld ( 33 ).

CONCLUSIONS
The location and dynamics of signaling lipids is paramount to their function, regardless of the precise signaling mechanism. We lack a clear understanding of both the exact cellular localization of many lipids and the driving forces that govern their movement. We ran atomistic simulations to determine the thermodynamics and kinetics of fl ip-fl op for three important lipids, as well as the free energies for the lipids' desorption in Raft and POPC bilayers. On the subnanometer length scale, we discovered some of the specifi c molecular determinants for the underlying free energy profi le, including hydrogen bonds, lipid membrane perturbations, and molecular geometry. Cholesterol, diacylglycerol, and ceramide have signifi cantly higher energy barriers for fl ip-fl op in the Raft bilayer than the POPC bilayer that translated into orders of magnitude across two membranes with diverse structures. Atomistic detail allowed us to determine both the number of hydrogen bonds to the lipid and its conformational fl exibility during translocation. The free energy barrier for fl ip-fl op was at the position where the number of intermolecular hydrogen bonds was zero for all three lipids. Our results support a solubility-diffusion mechanism for all three lipids in both bilayers, as opposed to a pore-mediated mechanism that has been observed for phospholipid fl ip-fl op ( 56 ). The three lipids are as long as half the bilayer thickness, so even at the bilayer center there are other contributions to the free energy, which argues the solubility-diffusion mechanism might be too simplistic. The Raft bilayer prevents conformational fl exibility, which is required during the process of fl ip-fl op and likely contributes to the higher energy barriers. These specifi c atomistic details would likely not be observed using continuum bilayer models, coarse-grained models, or most experimental techniques.

Lipid desorption
By comparing the free energies for lipid desorption between the Raft and POPC bilayers, we determined ⌬ G ex for all three lipids. Cholesterol had a large ⌬ G ex , which suggests a strong driving force for cholesterol incorporation into raft-like bilayer domains over nonraft (POPC) bilayer domains. Qualitatively, this matches the traditional view that suggests cholesterol preferentially interacts with saturated phospholipids and sphingomyelin. The preference of cholesterol for more ordered and rigid membranes with higher cholesterol concentration is important for passive cholesterol traffi cking from the endoplasmic reticulum with low cholesterol content, through the Golgi, to the plasma membrane ( 57 ).
In contrast to cholesterol, ceramide and diacylglycerol have relatively small ⌬ G ex , which indicates a similar preference for the POPC and Raft bilayers. The difference between ⌬ G ex for diacylglycerol and ceramide could be due to diacylglycerol having one unsaturated tail, whereas ceramide has two saturated tails. The low values of G ex suggest that at low concentrations ceramide and diacylglycerol would not have a strong preference to interact with raft domains over nonraft domains. It was shown that both ceramide and diacylglycerol are able to displace cholesterol from liquid-ordered domains and form more rigid domains ( 9 ). Our results suggest this behavior is concentration dependent, where a single ceramide does not displace a cholesterol from the domain, but rather a high concentration of ceramide alters the domain structure, forming a more rigid and ordered domain. It is possible cholesterol would have a lower preference for the ceramide domain and that ceramide would have a strong preference for the domain. Prior to signaling, ceramide and diacylglycerol molecules are at a low ambient concentration (approximately 0.1-1 mol%) ( 13 ) and would not cluster in rafts, which could be important for high-fi delity signaling. After enzyme action, such as sphingomyelinase, the local concentration of ceramide would increase, its chemical potential would likely change, and a signaling platform could form. slower translocation across the Raft bilayer. Slower fl ipfl op in Raft membrane can be understood from the stiffness of the bilayer, which prevents water from entering the hydrophobic interior to hydrogen bond with the translocating lipid and prohibits the conformational freedom of the lipid. Cholesterol was found to have a large preference ( ⌬ G ex ) for the Raft bilayer compared with the POPC bilayer, which suggests that cholesterol can be "trapped" in raft domains both laterally and transversely. Ceramide and diacylglycerol fl ip across the Raft bilayer more slowly than across the POPC bilayer, but they both have a relatively small ⌬ G ex between the two bilayers or nearly equal preference for both bilayers, which could be important for their signaling functions. Our simulations add to the physiochemical characterization of lipid membranes.