Federico
Elías-Wolff
ab,
Martin
Lindén‡
c,
Alexander P.
Lyubartsev
b and
Erik G.
Brandt
*b
aDepartment of Biochemistry and Biophysics, Stockholm University, Stockholm, Sweden
bDepartment of Materials and Environmental Chemistry, Stockholm University, Stockholm, Sweden. E-mail: erik.brandt@mmk.su.se
cDepartment of Cell and Molecular Biology, Uppsala University, Uppsala, Sweden
First published on 2nd January 2019
Cardiolipin is a non-bilayer phospholipid with a unique dimeric structure. It localizes to negative curvature regions in bacteria and is believed to stabilize respiratory chain complexes in the highly curved mitochondrial membrane. Cardiolipin's localization mechanism remains unresolved, because important aspects such as the structural basis and strength for lipid curvature preferences are difficult to determine, partly due to the lack of efficient simulation methods. Here, we report a computational approach to study curvature preferences of cardiolipin by simulated membrane buckling and quantitative modeling. We combine coarse-grained molecular dynamics with simulated buckling to determine the curvature preferences in three-component bilayer membranes with varying concentrations of cardiolipin, and extract curvature-dependent concentrations and lipid acyl chain order parameter profiles. Cardiolipin shows a strong preference for negative curvatures, with a highly asymmetric chain order parameter profile. The concentration profiles are consistent with an elastic model for lipid curvature sensing that relates lipid segregation to local curvature via the material constants of the bilayers. These computations constitute new steps to unravel the molecular mechanism by which cardiolipin senses curvature in lipid membranes, and the method can be generalized to other lipids and membrane components as well.
Cardiolipin (CL) is a lipid that constitutes 10–20% of mitochondrial membrane phospholipids,28 and is a major constituent in many bacterial membranes,29e.g., about 5% of E. coli phospholipids.30 CL deficiency is associated with numerous diseases,31–34 and CL is believed to stabilize respiratory chain complexes in mitochondrial membranes.35,36 CL is also involved in proton transport along the mitochondrial membrane surface.37 The unique structure of CL, which resembles two phospholipids joined together via the head groups, yields a conical shape that is the origin of its curvature sensing abilities. For example, CL localizes to the curved poles and septa of rod-shaped bacteria,26,38 and promotes polar localization of certain membrane proteins.39,40 Further, ATP synthase (which exploits energy released via transmembrane proton transport) associates with the highly curved edges of the mitochondrial membrane cristae,18 where the high curvature might also promote local enrichment of CL.
The majority of experimental studies of lipid and protein curvature sensing are based on fluorescence methods that monitor how sensor molecules partition between regions of different curvatures.8–11,16,20,24,25 There are several challenges with these experimental methods. One is the difficulty to ensure that the curvature sensor and the surrounding membrane are not influenced by dyes or other additives.41 For supported lipid bilayers,10,42 unwanted interactions with the solid surface must be avoided. Finally, fluorescence microscopy only provides the position and/or local density of the fluorescent labels,43 and therefore gives only indirect information about the structural basis of curvature sensing.
Fortunately, computer simulations have potential to provide additional details on the microscopic level, but face methodological and technological challenges of their own. Simulating curvature sensing and deformation in a flat bilayer patch relies on induced local membrane deformations which can be difficult to interpret.18,23,44–46 Curved membrane surfaces can be generated by simulating spheres47,48 or cylinders,49 but require equilibration of lipids and solute in different leaflets and compartments, which is only feasible for strongly coarse-grained membrane models. In this work, we use simulated membrane buckling50–60 and develop an elastic theory to overcome these difficulties. Our buckling method59 is based on a reaction coordinate (the arc-length) along the membrane profile, which presents a continuous range of curvatures in a single simulation, and is amenable for comparison to elastic membrane theories. In our simulations, the buckled membrane patch (Fig. 1a) consists of two equivalently deformed monolayers, which due to periodic boundary conditions does not divide the solutes into multiple closed compartments. Thus, the distributions of various species along the buckled membrane reflect their respective curvature preferences,53,54,56–60 and equilibration is limited only by in-plane lipid diffusion. We use the coarse-grained Martini model61–63 (Fig. 1b) to sample curvature dependent distributions and structural variations of CL mixed with the two other major E. coli lipids,64 1-palmitoyl-2-oleoyl phosphatidylethanolamine (POPE), and 1-palmitoyl-2-oleoyl phosphatidylglycerol (POPG). We find a strong curvature dependence on the lipid structure. The curvature dependent distributions are well explained by different intrinsic preferences for positive, zero, and negative curvature by POPG, POPE, and CL, respectively. The strength of these curvature effects is relevant to membrane structures such as mitochondrial cristae,18 and inclusion bodies resulting by over-expression of bacterial membrane proteins.12
A buckled bilayer is produced by assembling a flat bilayer in the xy-plane with a specified Lx/Ly ratio (Lx > Ly). After the bilayer is equilibrated, the bilayer is compressed along the x-direction, keeping the length of the simulation box in the y-direction fixed. The compression is performed by rescaling the simulation box instantly by a certain compression factor Lx/L = 1 − γ (where γ = 0 indicates a flat bilayer). The curved bilayer is energy minimized, re-equilibrated, and then simulated in production runs of 20–30 μs. We computed lipid autocorrelations to estimate sampling times, and found that single lipids are correlated for 120 ns in a pure POPE bilayer (details are given in Section S1 and Fig. S1 in the ESI†). Since the autocorrelation time is closely linked to lipid diffusion, we expect this result to approximately hold for all the simulated lipid systems. The first 1 μs was then discarded as equilibration in all simulations.
The simulated bilayers contained 512 lipid molecules, but with varying lipid composition. The charged systems were neutralized with counterions, and the amount of water in the systems corresponds to fully hydrated bilayers.72 The ratio between the lateral box lengths was Lx/Ly ≈ 4.5. Two single-component systems with 512 POPG and 512 POPE were included. The other simulations were multicomponent systems, with 384 POPE lipids, 0, 8 or 24 CLs, and the rest being POPG lipids. We used a CL model developed for Martini68 that corresponds to di-[1-palmitoyl-2-oleoyl-phosphatidylglycerol], corresponding to the most common type of cardiolipin found in E. coli membranes.73 The phosphate head groups (Fig. 1b) of POPG and CL are charged (−1). We used a symmetric head charge distribution for CL, with −0.5 charge units on each phosphate group, conforming to recent experimental data.74–76 We simulated two extra replicas of the system with 24 CL molecules (CL12): one with added salt (at 0.15 M), and one with a higher degree of buckling (γ = 0.4 compared to the other simulations at γ = 0.3). Since the presence of divalent cations can decrease the molecular surface area77 so that CL undergoes a lamellar-to-hexagonal (La-HII) phase transition,78 we used only monovalent ions (Na+ and Cl−) to neutralize and/or change the salt concentrations in the simulated systems. Table 1 summarizes the simulations carried out in this work. Note that all reported times are actual simulation times and have not been scaled with the semi-universal factor of 4 sometimes employed for the Martini model.79
POPG | POPE | CL | Na+ | Cl− | γ | L x ·Ly [nm2] | L z [nm] | Sim. time [μs] | |
---|---|---|---|---|---|---|---|---|---|
allPG | 512 | 0 | 0 | 512 | 0 | 0.3 | 113.4 | 18.7 | 20 |
allPE | 0 | 512 | 0 | 0 | 0 | 0.3 | 114.0 | 30.0 | 20 |
CL0 | 128 | 384 | 0 | 128 | 0 | 0.3 | 115.2 | 27.7 | 30 |
CL4 | 120 | 384 | 8 | 128 | 0 | 0.3 | 115.3 | 17.2 | 30 |
CL12 | 104 | 384 | 24 | 128 | 0 | 0.3 | 120.3 | 28.4 | 30 |
CL12s | 104 | 384 | 24 | 335 | 207 | 0.3 | 119.7 | 25.1 | 30 |
CL12b | 104 | 384 | 24 | 128 | 0 | 0.4 | 105.4 | 32.1 | 30 |
The membrane shape obtained with this buckling procedure is the two-dimensional analogue to the compressed elastic rod known as the Euler elastica.51,52 A theoretical description of the bilayer membrane as a thin elastic sheet leads, via a Helfrich functional model,80 to a parameterization of the elastica in terms of the arc-length parameter 0 ≤ s ≤ 1 (Fig. 1(a)). The free energy is56
(1) |
As an alternative, we have developed a numerical method which avoids special functions and efficiently extracts accurate values for X(s) and Z(s) from eqn (1).59 This method relies on Fourier expansions of X(s) and Z(s), with Fourier coefficients stored in lookup-tables for fast evaluation. The heart of the algorithm is to minimize the sum-of-square distances from the innermost lipid tail-beads (xj,zj) to their corresponding positions along buckled shape X(sj,γ), Z(sj,γ). The parameters in the minimization procedure are the offsets in the x and z directions, the sj (j = 1,…,N where N is the number of innermost tail-beads) which are the projected positions of the innermost tail-beads along the buckle in terms of s, and the compression factor γ (to prevent area fluctuations from biasing the fit). In practice, γ-fluctuations are only a few percent.56 Implementation details and speed and accuracy tests of this method can be found in ref. 59. An open source Matlab implementation of the Euler elastica is available within the mxdrfile package.81
The local curvature at position s is directly evaluated by the standard formula
(2) |
Fig. 1 shows a snapshot of a buckled bilayer, annotated with the fitted shape along the bilayer midplane. Position data for membrane components relative to the buckled profile can easily be extracted from an aligned trajectory. The projected s-positions of the innermost lipid tail-beads serve as the reaction coordinate to measure the curvature preference of each lipid component.
Two recent simulation studies58,60 of lipid curvature sensing in buckled membranes used a different analysis approach. Those studies employed an analysis method based on spline-fitting the lipid head beads and then extracting the local curvature by differentiation (i.e., by eqn (2)). Our analysis of the average bead densities, shown in Fig. 2a and Fig. S2 (ESI†), strongly suggests that the distance between the planes defined by the head groups vs. the bilayer midplane is independent of curvature, so that translating between these planes is only a matter of geometry. While the spline-fitting analysis is more flexible in terms of bilayer shape, our approach makes explicit use of elasticity to parameterize possible bilayer shapes. Further, since we have a theoretical reference for the buckled shape, we can align many frames to compute averages such as Fig. 2. More importantly, the well-defined particle positions sj along the buckled shape enables direct comparison to physical models developed using the tools of statistical physics.56,59 Thus, with our approach we can move beyond a descriptive analysis and directly validate or refute different model hypotheses. Since the reference shape of our membrane buckle is based on a model of solid physical merit (the Euler elastica52), our results furnish an enduring starting platform for future refinements.
Fig. 2 Heat maps of head- and tail bead densities in the xz-plane. (a) All lipid components for the CL12 simulation (Table 1). The red curve corresponds to the average fit (X(s; 〈γ〉), Z(s; 〈γ〉)). The blue curves are parallel to the first one, but displaced by 2.1 nm. (b) Same data as in panel (a) but for the CL component alone. Head bead densities have been weighted by factors 4 (POPG, POPE) and 8 (CL) to compensate for the number of tails per lipid, and the fact that the midplane contains tail beads from both monolayers. |
Geometrical effects of bilayer curvature with respect to lipid packing are apparent in Fig. 2a. While the bilayer innermost tail-bead density is almost constant, the head-bead regions corresponding to negative curvature are enriched, while the opposite is true for positive curvature. This effect is much more pronounced for cardiolipin as shown in Fig. 2b.
However, our simulations show asymmetric curvature dependence of the lipid distributions. This asymmetry is observed in Fig. 3a and c, which show the monolayer lipid density ρ(s) defined by the innermost tail-beads. The density function ρ(s), takes the perspective of the upper monolayer, with s = 0.5 corresponding to the region of maximal positive curvature, and s = {0, 1} corresponding to largest negative curvature. Contributions to ρ(s) from lower monolayer lipids are constructed by shifting their s-values by half a period. Fig. 3a and c show that the depletion around the region of minimal curvature is larger than the enrichment at maximal curvature. Also, the profile ρ(s) is flat around maximum curvature (s = 0.5). The bilayer midplane density (s) is depicted in Fig. 3d. This alternate function corresponds to the distribution of the innermost tail-beads of all lipids, without regard to which monolayer they occupy. Fig. 3d shows that this density does not deviate more than a few percent from (s) = 1, yet the deviations are systematic and consistent with the monolayer asymmetry. (s) is maximal around s = {0.25, 0.75}, where the bilayer is flat and lipid packing is tightest (densities are normalized as probability densities). Interestingly, we have previously seen the opposite trend for an even more coarse-grained three-bead lipid model.59
Fig. 3 The lipid distribution along the curved bilayers. (a) The monolayer density ρ(s) for all lipids as a function of s. (b) Curvature, K(s), (defined in eqn (2)) as function of s for the two γ-values used in this work. (c) Monolayer density ρ(s) for all lipids as a function of curvature K(s). Symbols correspond to different lipid compositions, degrees of buckling and salt concentrations (Table 1). The s-values correspond to the upper monolayer. The s-positions of the lower monolayer have been translated by half a period, according to ρ(s) = (ρupper(s) + ρlower(s + 0.5))/2. (d) Bilayer midplane density (s), which includes both monolayers. The sizes of the standard errors of the mean (SEMs) is smaller than the symbols in panels (a) and (c). In panel (d), SEMs are slightly larger than the symbols but are omitted for clarity. Fig. S5 in the ESI† shows a scaled-up version of this figures with SEMs. |
Fig. 4 Average order parameter 〈P2〉 as function of K(s) for saturated (open symbols) and unsaturated (solid symbols) tails. Panels correspond to (a) POPG, (b) POPE, and (c) CL, and symbols to lipid bilayers with different compositions, degrees of buckling, and salt concentrations (Table 1). Shaded areas show the standard errors of the mean (SEMs), as described in Section IIC. SEMs are only shown for CL0, CL12 and CL12b for clarity. SEMs are only noticeable in the regions of largest negative curvature in (a), where POPG is depleted, and are smaller than the symbols in (b). |
The presence of CL leads to a small increase in 〈P2〉 for the remaining lipids. We believe that this is due to the tails of CLs filling the larger available volume in regions of negative curvature. This mechanism is similar to how cholesterol can “plug the holes” in fluid bilayers that are rich in unsaturated lipids.84 Analogously, the ability for CL to stabilize protein–lipid domains has been ascribed to its ability to fill gaps in protein interfaces.85
Fig. 5 shows the normalized relative densities Φj of each lipid component j as function of local midplane curvature K(s). These densities are calculated as follows: From an aligned trajectory, the s-positions of the innermost tail-beads are collected to form the monolayer distributions ρj(s) for each lipid species j (Fig. S2 in the ESI†). ρj(s) is normalized to , where fj is the fraction of lipids of species j in the monolayer. The relative densities of each lipid component are computed by dividing ρj(s) by the corresponding all-lipid density ρ(s) (shown in Fig. 3a), i.e.,
(3) |
Fig. 5 Normalized relative densities Φj(s) as functions of curvature K(s) for each lipid component. Φj(s) is calculated as the density ρj(s) (Fig. S2, ESI†), divided by the corresponding all-component density ρ(s) (Fig. 3), and normalized as a probability density. The panels correspond to (a) POPG, (b) POPE, and (c) CL, and symbols correspond to the simulated systems as reported in Table 1. The shaded areas show the standard errors of the mean (SEMs), as described in Section IIC. SEMs are only shown for CL0, CL12 and CL12b for clarity, and are smaller than the symbols for all values of K(s) in panels (a) and (c). |
CL localization to negative curvature depends not only on its shape,60 but can also be affected by its headgroup charge.67,86 However, the ionization state of the CL headgroup has been heatedly debated in the literature (see ref. 87 and the references therein), partly because of its potential role in the molecular pathways for Barth syndrome33 and/or as a proton reservoir for proton-pumping respiratory enzymes.37 Recent FT-IR measurements,76 titration experiments,75 and electrokinetic/spectroscopic methods74 lend convincing support to CL being fully ionized at physiological conditions, with a symmetric charge distribution (−1 on each phosphate head group). We used a reduced-charge model of CL with −0.5 charge units on each phosphate group, to partially balance the overestimated electrostatic interactions in the Martini model (which uses relative permittivity εr = 2.5 with the polarizable water model).
We tested if increased charge screening can enhance CL curvature sensing in a simulation with 0.15 M salt concentration (CL12s), which is equivalent to that in blood plasma at physiological conditions. We found only minor changes to all calculated quantities (lipid distributions, chain order parameters, and enhancement ratios, as shown in Fig. 3–6) compared to their salt-free counterparts. The differences are most pronounced, but within 1–2%, at high positive and negative curvatures, and negligible at intermediate curvatures. These observations are consistent with earlier simulations58 that also showed a small effect on the degree of CL localization when the head group charge changed from −1 e to −2 e.
Fig. 6 Enhancement ratios for the two- and three-component simulations. ϕ+j(s) is the local relative density of lipid species j in regions of positive curvature (an innermost lipid tail bead with s ∈ (0.25, 0.75) corresponds to positive curvature if it resides in the upper monolayer), and ϕ−j(s) denotes the relative density of lipids in negative curvature regions. Black lines are least-squares fits to eqn (10). Shaded areas show the standard errors of the mean (SEMs), as described in Section IIC, for the CL0, CL12 and CL12b systems. |
Chen et al. have shown that 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine (DOPE) is more prone than CL to assemble into highly negative curvature structures in the Hexagonal II (HII) phase.88 Furthermore, CL curvature sensing depended on the presence of divalent calcium ions. This is most likely a characteristic of the HII phase, where CL curvature sensing is strongly suppressed by electrostatic repulsion between headgroup charges. The dominating factor in our simulations of the buckled lamellar phase is lipid packing frustration.
Following ref. 3 and 47, we model the curvature dependent energy per lipid molecule of species j in a monolayer with the harmonic expression
(4) |
S = −kB(logϕj(s) − 1), | (5) |
(6) |
(7) |
(8) |
(9) |
(10) |
A direct comparison of our simulations with experimental elastic constants is not possible, since we only extract combinations of products of the bending moduli and intrinsic curvatures for the lipid species. However, it is possible to do a consistency check on the product of elastic constants by assuming typical literature values for the bending moduli89 (25kBT) and area per lipids90 (0.7 nm2). With these numbers for all lipid species, the extracted products of bending modulus and intrinsic curvatures can be solved to yield the intrinsic curvature differences KPG − KPE ≈ 0.06 nm−1, KPG − KCL ≈ 0.26 nm−1, and KPE − KCL ≈ 0.17 nm−1. These differences are similar in magnitude to the experimentally estimated intrinsic curvatures of CL (−1 nm−1 to −0.2 nm−1)91 as well as other lipids.24,92
Finally, we emphasize that although Fig. 6 shows enhancement ratios versus curvature, the comparison is based on number densities along the buckled profile. Our analysis therefore relies on our ability to extract lipid arc-length positions sj from the simulations, which is also the reason for the unequal spacing of points along the K(s)-axes in Fig. 6.56,59
The shapes of the simulated buckled bilayers were in excellent agreement with predictions from elastic theory,51,52 and we observed no curvature dependent bilayer thinning for any of the lipid mixtures. We found that the curvature dependence of the lipid packing was asymmetric, with the density depletion being largest at the minimal curvatures. Our analysis shows that CL lipids strongly prefer membrane regions with negative curvatures. POPE lipids also prefer negative curvatures, but more weakly, and are outcompeted for the most negative curvatures by CL. POPG lipids show opposite behavior and are enriched in the regions on positive curvature. The relative preferences of CL and PE for negative curvature agree qualitatively with the observations of Boyd et al.58 They studied mixtures of POPE, POPC, and CL in a buckled geometry using the Martini model and a simulation setup similar to the one used here, but with a different analysis method.
Since lipid redistribution can stabilize curved membranes,93 we calculated the chain order parameters of the lipid tails along the curved surface to assess the curvature dependence of the lipid structures. We found the largest decrease in the chain order parameter in regions of maximal negative curvature for CL, where it is most concentrated. This is consistent with a geometric mechanism for curvature sensing, with cardiolipin “repairing” local membrane deformations in the curved structure.
Finally, we developed and validated a Helfrich-type elastic model that quantifies the competition among lipid components for the most favorable curvature region. The theory provides a simple relation between enhancement ratios47 for multicomponent membranes and the local curvature, and can be used to determine some of the membrane's material constants.
CL has theoretically been predicted to form clusters that localize to the poles of E. coli membranes once a critical concentration threshold is reached, partly driven by lipid–lipid interactions.27 We see no significant concentration dependence of our CL enhancement ratios (Fig. 6b and c), as expected for a model without lipid–lipid interactions (eqn (6)). However, while the nominal CL concentration in our simulations reach well above the predicted threshold of 1%,85 the absolute number of CL molecules may be too small to form clusters. We also probe curvatures more than an order of magnitude larger than those of bacterial cell poles, so it is perhaps not surprising that the curvature mismatch energy strongly dominates in our results. However, this regime is still of biological relevance for the extreme curvatures present for example in the mitochondrial inner membrane.18
CL is the signature lipid of mitochondria. Several proteins bind cardiolipin with high affinity, including respiratory complexes I, III, IV and V.94 The efficiency of ATP synthase is linked to high membrane curvature in mitochondrial cristae,95,96 and specific interactions between CL and the rotary mechanism of ATP synthase.97 CL segregation in the highly curved mitochondrial inner membrane may also enhance ATP synthase function by acting as a proton sink.98 The computations and analysis reported in this work are stepping stones for the theoretical explanation of curvature sensing by cardiolipin in lipid membranes at a molecular level, and can be generalized to study the general interplay between shape and composition in biological membranes.
Footnotes |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8sm02133c |
‡ Present address: Scania CV AB, Södertälje, Sweden. |
This journal is © The Royal Society of Chemistry 2019 |