The Molecular Structure of the Liquid Ordered Phase of Lipid Bilayers

Building the initial configuration

Two separate bilayer systems were built with the CHARMM-GUI, one with an L_o composition and one with an L_d composition.²⁹ The compositions were taken from published experimental results for a 1:1 ratio of DOPC:DPPC + 20 mol % chol at 298 K.¹³ The L_o system comprised 76 DOPC, 280 DPPC, and 156 CHOL, and the L_d system 240 DOPC, 116 DPPC, and 44 CHOL. Each bilayer was simulated in the constant particle number, pressure, and temperature ensemble (NPT) until the membrane area equilibrated. The initial configuration for the L_o/L_d system was then created by embedding a 5 nm diameter, roughly circular patch of the equilibrated L_o bilayer into the equilibrated L_d bilayer, and eliminating any overlapping lipids from the L_d bilayer (Figure S1). Care was taken both to maintain the correct composition of each region and to obtain a symmetric bilayer during this process. After embedding the L_o domain and symmetrizing the bilayer, the final number of lipids of each type was 174 DOPC, 114 DPPC, and 60 CHOL. This gives a bulk composition that is closer to 1:1.5 +20% DOPC:DPPC + CHOL, but since the miscibility transition temperature and compositions of each phase are only weakly dependent on the DOPC fraction,¹³ we expect this deviation to have a minor effect. A water box containing 16,261 waters and 0.150 M NaCl was added. The resulting system was allowed to relax under NPT conditions again, while weakly restraining the lipid headgroups to maintain the monolayers and also prevent water from entering the gap between the L_o and L_d regions. After the system area had relaxed, the restraints were removed, and system was relaxed for just over 25 nanoseconds of unrestrained NPT simulation. The end result of this procedure is shown in Figure S1, panel C. An initial configuration for a control simulation for the high temperature, mixed phase was built with the CHARMM-GUI, this time simply choosing the positions of all three components randomly.

Production simulation on Anton

Anton is a special purpose machine designed for high performance molecular dynamics calculations. Details on the Anton machine can be found elsewhere ³⁰. The equations of motion were integrated with the Verlet algorithm with a timestep of 2.0 fsec. A constant temperature and a pressure of 1 atm were maintained by the Martyna-Tobias-Klein method,³¹ with the pressure coupling effected every 240 fsec and the temperature coupling every 24 fsec. The temperature of the system initiated from a phase separated state was maintained at 298 K, the temperature of the homogeneous system was maintained at 328 K. Lennard-Jones interactions were truncated at 10.14 Å by a hard cutoff with no shift. Long range electrostatics were computed by the k-space Gaussian split Ewald method³² on a 64×64×64 point grid, with the parameters of the Gaussian chosen to yield a root mean squared error in the electrostatic force calculation of 0.18%. The duration of the T = 298 K system was 9.4 μsec, the duration of the T = 328 K system was 6.9 μsec, and the duration of the L_o control was 5.0 μsec.

Local composition analysis

We calculated the locally averaged area density of DPPC on a 20 × 20 grid (roughly 5 Å between grid points), counting DPPC phosphates within 10 Å of each grid point for each simulation snapshot. This quantity we call ρ(r), where r labels the grid point. To obtain a range for this observable that crosses zero at the average and has a range of order 1, we rescale this quantity by subtracting and dividing by the bulk average of ρ(r): ρ̃(r) ≡ (ρ(r)− 〈ρ〉)/〈ρ〉, where 〈ρ〉 is area density of DPPC averaged over the entire simulation cell. For this composition, 〈ρ〉 ≅ 0.57nm⁻². ρ̃(r) ranges from −1 (no DPPC) to a (positive) maximum value that corresponds to the highest observed local density of DPPC, roughly 4 for both temperatures.

Hidden Markov Model technical details

The HMM consists of two hidden states (putatively L_o and L_d) and 28 emission signals (the 28 different local lipid compositions). For each lipid in each frame of simulation (sampled every 0.239 ns), the local composition is measured. The time-ordered local composition is the emission signal for that unknown state. Between each frame, the lipid may change state (e.g., L_o to L_d) with a probability that does not depend on its history. There are four such so-called transition probabilities; each of the two states has a probability of changing or not changing state.

There are a total of 60 HMM parameters: 28 emission probabilities for each state and two transition probabilities for each state. The probabilities are constrained to sum to one, yielding 56 unique values to determine. These are determined by the Baum-Welch algorithm.^33,34 For some initial “guess” parameters, the algorithm computes the probability of observing the composition sequence given all possible sequences of states (summation over the astronomical number of paths is simplified by taking advantage of the Markovian character of the system via a forward-backward algorithm). The total expected number of emissions and transitions that occurred (summing over each possible path) are normalized and used as probabilities for the next iteration. Although an iteration is guaranteed to produce a more likely HMM, convergence to a global maximum is not guaranteed. Convergence of the parameters was checked by starting from ten initial guesses and selecting the best solution. In most cases a majority of the runs converged to nearly the same (most probable) model.

Once the HMM parameters are determined, the most likely state sequence is determined by the Viterbi algorithm,³⁵ which iterates through a sequence, saving (at each frame) the most likely state path to that frame from the start. Again, Markovian character is crucial to calculate the probability of a longer state sequence from the previous one.

A comparison of the HMM state assignment to the local density of DPPC and the thickness in each leaflet is shown for several configurations for T < T_m (Fig. S2) and for T > T_m (Fig. S3).

Calculation of deuterium order parameters and quadrupolar spectra

Deuterium order parameters were calculated at each position along the aliphatic chains of DPPC according to

where θ is the angle between the carbon – hydrogen bond vector and the membrane normal, and the average is taken over snapshots separated by 2.4 nsec. The lipids are divided into two populations for the calculation based on the Hidden Markov Model analysis. The left and right panels of Figure S4 present the order parameters for the T < T_m and T > T_m system, respectively.

Calculation of ²H NMR spectra

The simulated spectra were calculated according to the formula

where Δν_q(n) are the quadrupolar splittings for each methylene- and methyl group of sn-1 and sn-2 hydrocarbon chains, e²qQ/h=167 kHz, and S(n) are the order parameters from the simulation. A random distribution of orientations of bilayer normals to the magnetic field was assumed.

Extraction of S_CD values from experimentally measured ²H spectra

The order parameters of hydrocarbon chains of DPPC in the spectrum of coexisting L_d and L_o phases were extracted by choosing a pure L_d phase spectrum recorded at higher temperature that perfectly matched the L_d component in the L_d/L_o spectrum. After intensity adjustment, this L_d spectrum was subtracted to yield the spectrum of the pure L_o component. Order parameters were extracted from the L_d and L_o spectra by integration as described earlier.³⁶ A few resolved quadrupolar splittings were assigned directly. The quality of order parameter extraction was verified by comparing measured and simulated spectra of L_d/L_o phase coexistence.

Calculation of lipid-lipid pair contact frequencies

Lipid-lipid and lipid-cholesterol boundaries are mapped by identifying pairs that share a boundary in a Voronoi tessellation of a configuration. The sample of configurations then provides information on the statistics of pairwise contacts. The extent to which the observed interaction statistics deviate from a random mixture is found by comparison to randomized configurations, under the constraint that the composition remains fixed. Two null hypotheses were considered: Randomization over the entire system, and randomization within each phase.

Lipids mix laterally on 10 μsec timescales

An initial configuration for a L_o nanodomain was prepared by embedding a 5 nm diameter L_o region in a 10 nm × 10 nm L_d bilayer (Fig. S1). The compositions of the two phases were set to the experimentally determined compositions for a ternary mixture of DOPC/DPPC/cholesterol just below the miscibility transition.¹³ (Details on the system set up and simulation protocol are found in Supplementary Materials (SM).) Fig. 1 illustrates the extent to which the lipids initially in the L_o domain have mixed by the end of the simulation. This may be quantified by asking the question: what is the probability of finding an original L_o lipid at a distance r from another original L_o lipid? Comparing cumulative distribution functions for these lipids (CDF) averaged over the first 100 nsec of the simulation and the last 100 nsec, it is clear that the initial domain of L_o lipids has dispersed (Fig. S2). Furthermore, comparison to the CDF of an equal number of randomly selected lipids shows that by the end of the simulation the initially L_o lipids have mixed (dotted line in Fig. S2), indicating that the simulation time-scale is sufficient to sample compositional degrees of freedom.

Local composition reveals a stable domain for T < T_m

Experimentally, it is known that the L_o and L_d regions differ in composition, with the L_o region comparatively enriched in DPPC and cholesterol. Whether a particular lipid belongs to the L_o or the L_d region should therefore be encoded in the composition of its local neighborhood. A neighborhood enriched in DPPC and cholesterol is more likely to be in the L_o phase than L_d, but how much more likely?

This was answered by developing a hidden Markov model (HMM) for the local membrane phase. The local composition in the vicinity of each lipid (determined by the six nearest lipids in the same leaflet, including itself) is recorded in a time series (one for each lipid). This “observable” data is used to define two “hidden” states (putatively L_o or L_d,). The parameters of the model are the probabilities of each hidden state to have a certain lipid composition, and the probabilities for a hidden state to change to the other in one time step. The values of the unknown probabilities were determined by the requirement that they maximize the likelihood of the observed sequence of 39,272 compositions for each of the 348 lipids using standard methods for solving HMMs. The HMM can also be used to define the hidden state at any point in the bilayer, given a local lipid composition. (See Methods for complete details.)

The result of the HMM analysis is a map of each leaflet of the bilayer onto L_o and L_d phases; a representative configuration is shown in the top panels of Fig. 2. The bilayer phase determined by the HMM analysis is consistent with a direct analysis of the local composition and the bilayer thickness (middle and bottom panels of Fig. 2, and S3), both of which are markers of the L_o phase. It is clear that despite extensive mixing of the lipids, a local region with features consistent with the L_o phase persists throughout the simulation.

The simulation admits direct observation of the dynamics of nanoscale L_o/L_d coexistence, included as a web enhanced object, with the rendering explained in the caption to Figure 4. (“lowT_small.mpeg”) The persistence of the L_o region is apparent, as is the exchange of lipids between the regions. Transient hexagonal structures are observed within the L_o region, discussed in more detail below.

In contrast, a simulation of the homogeneous phase above the miscibility transition temperature reveals compositional fluctuations that at times appear L_o-like, yet occur on a much faster timescale (Fig. S4 and movie available as a web enhanced object “highT_small.mpeg”). Despite a similar structure in the lateral composition, the T > T_m data have distinctly different structure in the hydrocarbon chains, described next.

Deuterium quadrupolar splittings confirm L_o/L_d coexistence

Based on the assignment of each lipid into either the L_o or L_d phase, we computed the C-D order parameters of DPPC hydrocarbon chains (S_CD) at each position for both temperatures (Fig. S5). S_CD reports the average alignment of the carbon-hydrogen covalent bonds relative to the bilayer normal, and is proportional to the quadropolar splitting of resonances observed in a ²H NMR experiment on samples containing DPPC with perdeuterated hydrocarbon chains. The quadrupolar splittings therefore distinguish the L_o and L_d phases by virtue of the more ordered chains in the L_o phase—indeed, this feature is the origin of the word “ordered” in “liquid ordered phase.”

In order to facilitate comparison to the experimental data, ²H NMR spectra were computed from the simulated order parameters as described in Methods. The simulated and measured spectra are shown in Fig. 3. The central resonances corresponding to the terminal methyls (red arrows in the bottom panel of Fig. 3) show a characteristic triple peak structure below T_m, but yield a single splitting above T_m. In the L_o phase (T < T_m), the lipids are oriented parallel to the bilayer normal with more ordered chains. Consequently, the sn-1 chain protrudes more deeply into the bi-layer center than the sn-2 chain, yielding distinct splittings and producing the two outer pairs of peaks below T_m. In the L_d phase the chains are disordered and the terminal methyls produce a single splitting corresponding to the innermost pair of peaks. The resonances arising from the rest of the carbons below T_m overlap and are difficult to distinguish on a carbon-by-carbon basis, but they yield two distinct populations, indicated by the two broad shoulders in the wings of the spectra. By comparison, the high T data in Fig. 3 are consistent with a single population of lipids without any components that are distinguishable on the NMR timescale of about 10 μsec. Table 1 compares the average order parameters from simulation and experiment. The chains of the L_o phase are slightly less ordered on average than reported by the experiment, most likely due to the high fraction of lipids at the L_o/L_d boundary in the T < T_m simulation. Indeed, the average order parameters observed in a control simulation of the L_o phase are in excellent agreement with the experimentally measured values. The compositions of the two phases also agree with those measured experimentally ¹³ (Table 2), except for the balance of DOPC to DPPC in the L_d phase.

Substructure composed of aligned alkyl chains with locally hexagonal order is observed in the L_o phase

Visualization of the centers of mass of the hydrocarbon chains and cholesterol reveals a local, approximately hexagonal substructure within the L_o phase, comprised of regions of densely packed DPPC chains, shown in Fig. 4A and in the web enhanced object discussed above. Similar substructure is observed in the L_o control simulation (Figs. 4B, 4C). A movie corresponding to Fig. 4B is available as a web enhanced object (“LoControl_small.mpeg”). These regions are composed of lipids with aligned saturated chains oriented parallel to the bilayer normal, as shown in the molecular graphic (Fig. 4C). They have sn-1- and sn-2 hydrocarbon chains with distinctly different quadropolar splittings of terminal methyl groups as measured for the L_o phase. The substructure is nonetheless fluid — the order does not extend over long range, it is dynamic, and analysis of the lifetimes of saturated chain contacts shows that contacts deepest in the hexagonal regions of the L_o domain have a lifetime of approximately 400 nsec (Fig. 5). Despite the compositional fluctuations that mimic the L_o phase, aligned chains, persistent contacts, and locally hexagonal order are not observed above T_m.

Because the simulation timescale admits thorough sampling of the compositional degrees of freedom, the statistics of lipid interactions may be addressed directly; here the frequencies of lipid-lipid contacts are analyzed, as measured by boundary lengths in a Voronoi tesselation. These boundary lengths are compared to the expected lengths for random distributions (ideal mixing) over the entire system and within each phase (Table 3, details in Methods). In every case, nearest neighbor cholesterol-cholesterol interactions are strongly disfavored, consistent with previous simulation results for binary mixtures.²⁰ Within the L_o phase, the local hexagonal structure is reflected in the frequency of DPPC-DPPC contacts, which are enhanced relative to ideal mixing. Cholesterol boundaries were distingiushed between the smooth and rough faces (S and R in Table 3). Within the L_o region, cholesterol-DPPC contacts favor the smooth face, while cholesterol-DOPC contacts favor the rough face. Taken together, these data indicate that cholesterol mediates the boundary between the locally hexagonal regions in the L_o phase and the interstitial, more disordered regions, and that it does so by virtue of its two (α and β) faces.

²H NMR quadropolar splittings

Distinct splittings of the sn-1 and sn-2 DPPC terminal methyls are produced by very ordered hydrocarbon chains aligned parallel to the bilayer normal. This is a distinguishing characteristic for detection of L_o phases. Although the hydrocarbon chains have high order, characteristic of chains with few if any gauche defects, the spectra still indicate rapid reorientation about the bilayer normal of lipids as individual molecules or small clusters. Furthermore, for resolved resonances only one set of quadrupolar splittings per methyl- or methylene group is detected, indicating that the chains explore their entire configurational space on the NMR timescale (correlation times less than 10 μsec). This imposes strong limits on aggregate size.¹³ In summary, the experimental spectra are in excellent agreement with organization of DPPC in small, highly mobile clusters with hydrocarbon chains packed in a hexagonal lattice and rapid exchange of DPPC molecules.

¹H magic-angle spinning NMR

The chain methylene resonances of DPPC in the L_d phase have a typical linewidth of 50–100 Hz. They broaden by a factor of ten in the L_o phase compared to L_d. This is the result of greatly reduced trans-gauche isomerization of chains, consistent with the arrangement of DPPC in mostly trans hexagonally ordered regions.³⁷

NMR diffusion measurements

Mean squared displacements are reduced by roughly a factor of 3 in L_o compared to L_d.³⁸ This is consistent with the slower (but still liquid-like) diffusion observed in L_o, clearly evident in the movie of the L_o phase accompanying Fig. 4.

Small angle X-ray scattering

Bilayer thickness is increased and area per lipid decreased in L_o,³⁹ due to regions of packed DPPC chains.

Small angle neutron scattering

Nanoscale domains with hexagonal order have recently been reported in binary mixtures of DPPC with cholesterol.⁴⁰

The composition and cholesterol dependence of the L_o phase evoke nanoscale size raft regions in the cell membrane, and yet membrane components that partition to raft regions in cell membranes typically do not partition to the L_o phase in bilayers.^6,41 Since many putative raft components have a preference for solvation by cholesterol, within the L_o phase they should partition to the cholesterol rich regions found at the boundaries of the densely packed hexagonal DPPC clusters. In the macroscopic L_o phase observed below T_m for mixtures of DPPC, DOPC, and cholesterol, there is comparatively little cholesterol rich area available to solvate such components. In the resting cell membrane, however, raft regions are known to be quite small,^9,10 consistent instead with nanoscale clusters of ordered chains, solvated by a cholesterol rich boundary region. These nanoscale clusters may be driven by T > T_m composition fluctuations, stabilized by immobilized “defects”,⁴² or they may be a result of reduced line tension from replacing DOPC by POPC.¹⁴ Since many small nanoscale clusters have an increased total boundary length relative to a macroscopic phase, components favoring cholesterol would find more such regions accessible. Furthermore, the size and balance of substructure within such regions (and therefore the area available to solvate raft components) in the cell membrane may well be shifted by several mechanisms, including coupling to the cytoskeleton,^42–44 proximity to the miscibility transition,^11,13,45 or active maintenance by the cell.^46,47 Indeed, crosslinking of components that partition into the densely packed substructure provides a natural mechanism to trigger small clusters to coalesce into a macroscopic domain.⁴¹ Thus, understanding precisely how the L_o substructure is altered by such mechanisms presents itself as a promising line of inquiry to resolve the connection between lipid bilayer and cell membrane lateral structure.

But why then is cholesterol a necessary component for liquid-liquid coexistence? Beginning with early work on binary mixtures of cholesterol and saturated phospholipids,^48–50 it has been hypothesized that the structure and thermodynamics of the L_o phase are driven by direct interactions between cholesterol and saturated lipids. This idea has been elaborated by McConnell and Radhakrishnan in terms of “condensed complexes.”⁵¹ The condensed complex model departs from regular solution theory by including a fourth species to model a stoichiometric condensed complex of cholesterol and (in this case) DPPC. An alternate view that proposes a longer range, lattice-like substructure in the L_o phase of binary mixtures containing cholesterol has been proposed by Chong and coworkers, with the lattice composed of condensed complexes of lipid and cholesterol.⁵²

In contrast, the present results indicate that cholesterol disrupts what would otherwise be gel-fluid coexistence, not just by direct interactions with the lipids, but indirectly, because cholesterol prefers an organization that is incommensurate with the local order preferred by DPPC in the L_o phase. Consistent with this idea, Martinez-Seara, et al.²⁰ recently showed that when in binary mixtures with DSPC, cholesterol prefers to interact with other cholesterols in next-nearest neighbor solvation shells at three-fold symmetric orientations, but the acyl chains of DSPC lie at roughly five-fold symmetric positions around cholesterol. Consequently, the rotational symmetry of cholesterol in L_o regions appears to be only partially compatible with hexagonally ordered DPPC chains. One may then view the condensed complex model as a way to incorporate non-mean field density fluctuations into regular solution theory. The present data indicate that these non-mean field fluctuations are not stoichiometric complexes of cholesterol with DPPC, but, rather, DPPC with itself. Indeed, more than 20 years ago Sankaram and Thompson suggested, on the basis of ²H NMR spectra and electron spin resonance measurments of binary mixtures, that the data could be explained not by stoichiometric complexes but by time averaged cholesterol positions.^53,54 This hypothesis could be tested by calculations similar to those presented here, replacing cholesterol by other sterols, such as lanosterol. Published 2H NMR⁵⁵ and quasi-elastic neutron scattering⁵⁶ studies report sterol-dependent differences in the ordering and dynamics of the lipid chains.

In summary, the L_o phase of a mixture of cholesterol and two lipids is shown to be itself inhomogeneous. Lateral segregation within the L_o phase is observed, with regions of hexagonally packed saturated chains separated by insterstital regions enriched in cholesterol and unsaturated chains. The observed substructure explains existing experimental data, and provides a focus for future efforts aimed at understanding the molecular scale structure of cell membranes.