7194
J. Phys. Chem. B 2007, 111, 7194-7200
Polymorphic Phase Behavior of Cardiolipin Derivatives Studied by Coarse-Grained Molecular Dynamics Martin Dahlberg* DiVision of Physical Chemistry, Arrhenius Laboratory, Stockholm UniVersity, S-106 91 Stockholm, Sweden ReceiVed: March 10, 2007; In Final Form: April 23, 2007
Cardiolipin (CL) is a negatively charged four acyl chain lipid, associated with energy production in bacterial and mitochondrial membranes. Due to the shape of CL, negative curvatures of aggregates are favorable if the charges in the head group can be reduced. The phase polymorphism of CL, and of associated derivatives with 2, 3, 4, or 5 chains, has been determined previously and offers a model system in which micellar, lamellar, and inverse hexagonal phases can be observed. We present an extension to a previously established coarse-grained molecular dynamics model with the aim of reproducing the different CL phases with two adjustable parameters: the number of acyl chains and the effective head group charge. With molecular dynamics simulations of large lipid systems, we observed transitions between different phases on the nanosecond to microsecond time scale. Charge screening by high salt or low pH was successfully modeled by a reduction of phosphate charge, which led to the adoption of aggregates with more negative curvature. Although specific ion binding at the interface and other atomistic details are sacrificed in the coarse-grained model, we found that it captures general phase features over a large range of aggregate geometries.
Introduction Cardiolipin (CL, Figure 1) is a naturally occurring tetraacyl phospholipid, found in bacterial and mitochondrial membranes. In human mitochondria, the inner mitochondrial membrane is especially rich in CL (9.2% by mole),1 and because CL is anionic, the membrane is negatively charged. The importance of CL for mitochondrial function can be explained on several levels. CL is virtually the only negatively charged lipid in the inner mitochondrial membrane,2 where a chemical gradient of protons is used in the respiratory chain to restore the fundamental energy currency of the living cell, ATP. Indeed, CLproton interactions have been linked to charge carriage in oxidative phosphorylation.3 Additionally, CL coseparates with most respiratory chain proteins, the efficiency of which are dependent on CL levels.4-9 An example of interactions between CL and cytochrome c (cyt c) has been investigated recently by Gorbenko et al.,10 who determined quantitatively the electrostatic conditions under which CL and cyt c form a lipid-protein complex. On a higher level, problems with CL metabolism due to mutations which decrease the ability to maintain proper CL levels can lead to heart failure, myopathy, and growth retardationsclinical features of Barth syndrome.11-13 The two most distinguishing features of CL are the four acyl chains and the potentially large negative charge of the phosphatidyl head groups. Both of these properties have a strong influence on the phase preference of the lipid. CL can readily form either lamellar (LR) or inverted hexagonal (HII) structures, depending on the pH and salt concentration and valence, as has been shown in a number of studies.5,14-17 Powell and Marsh18 studied CL derivatives with varying acyl content, and determined the phase diagram for CL derivatives with two to five chains per head group at several salt concentrations. A sequence from * E-mail:
[email protected].
Figure 1. Atomistic and coarse-grained representation of cardiolipin (CL4).
a micellar phase at low acyl and salt content via lamellar at intermediate acyl content to inverted hexagonal at high acyl and salt content was found. These changes in the preferred phase can be understood in terms of the dimensionless packing parameter, which relates the geometry of the lipid head group to the hydrophobic tail length and volume.19,20 The packing parameter is a valuable tool for understanding the connection between lipid geometry and aggregate geometry. For explicit predictions, however, the packing parameter is less useful, in part because of the difficulty in finding reliable values of the head group area, which can change as a result of changes in hydration, ion screening, and phase geometry. Potentially, both phosphate groups of CL can carry negative charges. However, these phosphate groups have very different pKa values (pK1 < 4.0 and pK2 > 8.0), which has been attributed to an intramolecular hydrogen bond with the 2′-hydroxyl group found on the glycerol linker in between the two phospholipids.3 The degree of protonation in mixed lipid systems is affected not only by the pH of the solution but also by the pH of the concentration of CL10 together with the ionic strength and valence of the solution. It should also be kept in mind that the pKa values measured in homogeneous solutions are slightly different from those found at the aggregate surface due to differences in the
10.1021/jp071954f CCC: $37.00 © 2007 American Chemical Society Published on Web 06/02/2007
Polymorphic Phase Behavior of CL Derivatives electrostatic potential and hydration,21,22 factors which make precise determination of the protonation state more involved. Recent work with surface force measurements emphasized that factors other than that of the pH in the solution affect the charge density of CL.23 It is generally accepted that polymorphic properties of lipids play a central role as intermediate structures in membrane fusion, as well as in protein insertion and translation.14,24,25 Because of the strong link between the property of the individual lipids constituting the membrane and the polymorphism, it is important to study these processes on a molecular level. Molecular dynamics (MD) computer simulations of lipids have successfully provided details on the atomic and molecular level, and today complement, for example, NMR and X-ray data in the study of membrane structure and dynamics. Coarse-graining from the atomistic scale to mesoscopic scales has recently gained much popularity. The procedure is a step in the effort to reach larger system sizes and longer time scales, and achieves that goal by successively decreasing the structural degrees of freedom. Adjacent atoms are joined into single interaction sites, and atomic details and absolute dynamics are sacrificed for speedups of several orders of magnitude. The potentials can be Lennard-Jones interactions, which can be parametrized against experimental or atomistic simulation data. Alternatively, effective potentials from previous atomistic simulations or distribution functions from scattering data can be used. Notable examples are the inverse Monte Carlo method,26,27 the force-matching method,28,29 and the Boltzmann inversion method.30 Combinations of Lennard-Jones interactions and numerical potentials derived from atomistic simulations have also been used.31,32 Coarse-grained (CG) methods and some recent studies have been summarized in the reviews by Venturoli et al.33 and Nielsen et al.34 Following the work of Smit et al.,35 Marrink and co-workers have developed a general CG lipid force field.36,37 In their scheme, CG particles interact via pairwise Lennard-Jones potentials, and Coulomb interactions with a cutoff. One CG particle replaces four to six heavy atoms in all-atom models, and bonded interactions between CG particles (bonds and angles) are described by weak harmonic springs. Zwitterionic lipids were previously simulated with the resulting force field, and the phase behavior was shown to agree well with experimental results.37 Recently, the same model has been used to study vesicle fusion,38 lipoprotein behavior and assembly,39,40 and membrane protein behavior.41 CG simulations are thus useful tools for lipid systems. At a significantly higher computational cost, large-scale atomistic MD simulations are however still used to gain more detailed insight into membrane structure and dynamics.42-45 In the present study, MD computer simulations were applied to a set of CG models of CL derivatives: dilysocardiolipin (CL2), monolysocardiolipin (CL3), cardiolipin (CL4), and acylcardiolipin (CL5). We vary the number of acyl tails together with the pH to induce phase transitions between micellar, lamellar, and inverted hexagonal phases. Our work is largely inspired by the 31P NMR and X-ray diffraction studies by Powell and Marsh.18 In addition to the biological relevance of CL4 and CL3, these lipid derivatives give an opportunity to study physicochemical problems such as the connection between charge and aggregate geometry. Because direct modeling of pH in molecular simulations is difficult both due to the low concentration of protons and hydroxide ions at physiologically relevant pH and due to the lack of atomic resolution in the CG model, another approach was used. We assume that a change in pH is reflected in the
J. Phys. Chem. B, Vol. 111, No. 25, 2007 7195 protonation state of the two phosphate groups of each CL derivative. Charge-charge interactions are thus modified to reflect the pH. Although charge interactions are included in the zwitterionic lipids for which the CG force field have been optimized, the present study is, to the best of our knowledge, the first in which CG simulations are carried out with charged lipids. Because of the long-range nature of charge interactions, special attention needs to be paid to the cutoff distance. Electrostatic interactions between head groups and ions in cardiolipin-rich lipid phases are most likely important for the structure of the lipid phase, and the distribution of ions close to the interface. It is also likely that interactions such as these constitute one of the reasons for keeping a high cardiolipin concentration in the inner mitochondrial membrane, where proton gradients are coupled to ATP synthesis.3 The ion distribution close to and at the interface is easily calculated from simulations. Although atomistic effects such as ion binding in the head group region46-48 are not visible at the CG level of detail, the main features of ions close to a charged surface should be preserved. Such systems are generally described well by the Gouy-Chapman theory, in which the ion distribution of a charged surface with counterions decays smoothly as a function of distance from the surface. Methods General Model. The CG force field, based on the work of Marrink et al.,36 was used. CL derivatives (Figures 1 and 4) were built from two (DOPE) lipids, where both ethanolamine groups were replaced by one glycerol particle, linking the phosphatidyl groups. Chains were removed or added as needed to make CL2, CL3, and CL5. The head group geometry (PO4glycerol-PO4) was enforced with a harmonic angle bending term with a force constant of 25 kJ/mol/rad2 and a potential minimum at 180°, in order to model the essentially linear geometry found in the crystal structure of CL. Glycerol groups without ester links and charge neutralized PO4 groups were changed from the “Na” and “Qa” types, respectively, to the “Nda” particle type36 to reflect hydrogen bond donor and acceptor ability. Monounsaturated oleoyl chains were chosen for all lipids in order to be able to use the existing oleoyl chains in the force field. The unsaturated bond in the oleoyl chain was modeled with a bond angle of 120° between the central three particles in the chain. In order to ensure that the liquid crystalline phase was simulated, all simulations were carried out at 293 K, substantially higher than the melting temperatures measured for the related CG models of dioleoylphospatidylethanolamine and dioleoylphospatidylcholine (DOPE and DOPC).36 In order to examine the effects of the double phosphate head group, we also constructed a dioleoylphosphatidylglycerol (DOPG) model based on DOPE as a reference system. Charges. The partial charges of the two PO4 groups were set either as -0.7e or as 0.0e, in an attempt to model titration from neutral to acidic pH. Effective partial charges of the anionic head groups are a function of the pH in the system. Full (-1.4e) and half (-0.7e) charge are defined as “neutral pH” and “low pH”, respectively, and represent the two lipid charge states investigated here. For the case of half charge, it can be argued that at the time scale studied (ns to µs) the charges are evenly distributed over the two phosphate groups (i.e., -0.35:-0.35 rather than 0.0:-0.7). However, using full charges makes it easier to use the original force field. Charge neutrality was enforced by replacing water particles with sodium ions bearing a +0.7 charge (as defined in the original force field36). The LJ interaction between water and phosphate was increased from 5
7196 J. Phys. Chem. B, Vol. 111, No. 25, 2007
Dahlberg
TABLE 1: Summary of Simulation Conditions for the Main Simulations in This Work model CL2 CL3 CL4 CL5 DOPG
charge/lipid
effective pH
lipids
water/lipidb
initial structurec
-1.4 -1.4 -0.7 -1.4 -0.7 -1.4 -1.4 -0.7 -1.4 -1.4 -0.7
neutral neutral low neutral low neutral neutral low neutral neutral neutral
500 8000 1000 864 864 384 864 864 864 864 124
37 20 18 8 8 45.7 8 8 8 12 45.7
isotropic isotropic isotropic bilayer (6) bilayer (6) bilayer (1) bilayer (6) bilayer (6) bilayer (6) bilayer (6) bilayer (1)
final structurea M M L L L L L HII HII HII L
simulation time (ns)
comment
200 60 200 400 160 400 160/360 160/100 230 2000 320
unbranched micelles branched micelles cross-linked bilayers
disordered Figure 5
a M, micellar; L, lamellar; HII, inverse hexagonal. b CL counted as two lipids, ion particles counted as four water molecules. c Bilayer (X) corresponds to X bilayers in a lamellar stack.
to 5.8 kJ/mol, in accordance with the changes reported by Marrink and Mark.37 A 2.0 nm cutoff with a shift from the origin was used for the electrostatic interactions. Because long-range interactions are crucial for charged amphiphiles, different cutoffs were tried. The relative permittivity was set to 20, as in the original force field, and Lennard-Jones interactions were cut off at 1.2 nm. Simulation Details. Temperature, 293 K, and pressure, 1 bar, were controlled using a Berendsen49 thermostat and anisotropic barostat, both with a coupling parameter of 1 ps. Compressibility in all directions was set to 6 × 10-5 bar-1, corresponding to the compressibility of CG water.36 To investigate the hexagonal and micellar phases, fully flexible triclinic boxes were used in order to not bias the system toward the lamellar phase. A time step of 40 fs was used. The GROMACS (3.3 and 3.2.1) package50 was employed for preparation, simulation, and analysis. The time scale in the simulations was matched to the diffusion coefficient of water,36 by multiplying the time in simulation with a factor of 4. Composition. Lipid phase preference is not only dependent on the geometry (e.g., acyl chain content) and head group interactions but also dependent on the water content. The water contents used in our simulations are given in Table 1, where other parameters for the main simulations are also summarized. In order to adhere to the experimental protocol of Powell et al.,14 where a 14 wt % minimal water content for the CL4 HII phase was measured, we used a water to lipid ratio of 4 particles of water per CL4 (27 particles). This ratio corresponds to 15 wt % water in the case of CL4. To facilitate comparison with other membrane lipids, we give values of water molecules per lipid with cardiolipin counted as two lipids. Then, 15 wt % water corresponds to eight waters per lipid. In the calculation of water content, sodium particles were assumed to include four water molecules in the solvation shell, but it should be noted that, given the coarse-grained nature of the model, exact water contents are not available. Higher water contents were used for simulations from isotropic mixtures, in order to avoid semistable intermediates. Initial coordinates were generated from smaller self-assembled systems by stacking several boxes together. In the case of CL2, simulations were started from isotropic random mixtures. In order to assess the phase stability for CL3 and CL4, starting structures were built from simulations of CL4 at low pH. For CL3 simulations, one oleoyl chain per CL4 was removed and the resulting system was energy minimized. Results Treatment of Charges. In order to test the assumption of smoothly decaying ion distributions, bilayers with 96 CL4 lipids
were simulated at a water content of 45.7 waters per lipid. The high water content prevents strong interactions between adjacent bilayers. In Figure 2, the counterion density distributions for the CL4 bilayer are plotted as a function of distance from the bilayer center. The density distributions were calculated as the number density of the charged particles along the normal of the bilayers. Four different treatments of the Coulombic interactions were compared: particle mesh Ewald (PME) and cutoffs at 1.2, 2.0, and 2.4 nm, respectively. For the three most long-range treatments of electrostatics, the density distributions are of similar shape, and the intensity of the peak close to the membrane differs at most about 14%. With a cutoff of 1.2 nm, a significant shoulder appears in the distribution just outside the main peak (at (2.5 nm in Figure 2), reflecting an unphysical correlation close to the cutoff range. Since these model membranes bear a high surface charge, the Gouy-Chapman length is small, implying that a large part of the ions will be contained close to the membrane surface. However, because the charge distribution of the head groups is also quite wide along the normal to the membrane surface, the counterion distribution is widened. As the cutoff is increased, the unphysical shoulder in the distribution disappears. This is most likely because the bulk of the ions are contained within 1 and 2 nm of the membrane. The radial distribution function between PO4 groups was calculated (data not shown) to ensure that no spurious peaks appear at the cutoff distance. Indeed, no such peaks were found for any cutoff. A 2.0 nm cutoff was used in all further simulations. The 2.4 nm cutoff and the PME method were avoided for reasons of computational efficiency. The ion distributions for these two methods differed little from that of the 2.0 nm cutoff. By changing the force field, which was parametrized with a 1.2 nm cutoff for charges, the balance of forces is changed. In order to assess the impact this has on equilibrium bilayer structure, the area per lipid was calculated for the different electrostatics treatments in simulations lasting 400 ns each. The area per lipid was calculated from the box area in the plane parallel to the membrane surface, divided by the number of CL molecules per leaflet. In order to facilitate comparisons with DOPG and other two acyl chain lipids, the area for CL has been multiplied by 0.5. A difference of about 1% was found between the average areas per lipid, with the 2.0 nm cutoff at 0.653 ( 0.003 nm2/lipid and the 1.2 nm cutoff at 0.646 ( 0.002 nm2/ lipid. This difference is small and probably not crucial for the head group induced changes in phase preference. In order to assess charge screening without reducing PO4 charge directly, we carried out a 600 ns simulation of a CL4 bilayer with 364 lipids with approximately 1.5 M monovalent
Polymorphic Phase Behavior of CL Derivatives
Figure 2. Sodium ion number density as a function of distance from the bilayer center and range of electrostatic interactions. The density has been averaged over the two sides of the bilayer.
Figure 3. Order parameter, P2, of bonds in DOPG (circles) and CL4 (+). GL* denotes the linking glycerol particle in CL4 and the terminal glycerol in DOPG.
salt added. The area per lipid decreased only slightly (0.640 ( 0.003 nm2). Although a large proportion of the head groups should be screened by such salt concentrations,18 the area per lipid change was about 1%. It is likely that when enough salt is added to a real CL-water system, the additional dehydration18 of the head group region together with charge screening allows for more negative curvatures, and the transition into the hexagonal phase is possible. From the small change in area per lipid with 1.5 M added salt, we conclude that screening of the head groups is too low to allow for salt induced transitions with the CG model. The inability to accurately model salt effects on the lipid head groups is a shortcoming of the model, which was developed for use close to physiological conditions, with salt concentrations in the 0.1-0.2 M range. We should not expect the model to be valid at salt contents in the 1-4 M range and therefore focused the simulations in the present work on changes in the head group charge directly. Bilayer Properties of CL4 and DOPG. Lamellar patches of CL4 and DOPG were simulated to establish differences between the two closely related molecules. The area per lipid was slightly smaller for cardiolipin (0.643 ( 0.003 nm2) than for DOPG (0.658 ( 0.007 nm2), and was in between two experimentally determined values, 0.71 nm2 14 and 0.6 nm2.10 For a fully saturated tetramyristoyl-CL, an area per lipid of 0.5 nm2 was recently reported,51 which is reasonable given the lack of kinks in the acyl chains. The difference between DOPG and CL4 in the simulation can be explained from the fact that
J. Phys. Chem. B, Vol. 111, No. 25, 2007 7197 PO4 groups are constrained through the connecting glycerol group in cardiolipin, a constraint that two negatively charged groups in neighboring DOPG lipids do not have. DOPC is similar to DOPG but zwitterionic, and the coarse-grained model for DOPC has an area per lipid of 0.67 nm2 at 300 K.36 Charged lipids commonly have areas per lipid similar or below those of similar zwitterionic lipids, a phenomenon observed for, e.g., DMPS, DPPG, DOPS, and POPG.48,52 The second-rank order parameter in the lamellar CL4 and DOPG phases is defined by P2 ) 1/2〈3 cos2 φ - 1〉, where φ is the angle between the vectors joining consecutively bonded particles and the z-axis of the box. The brackets denote averaging over molecules and time. Order parameters for bond vectors in DOPG and CL4 were compared and showed a slight increase in the chain order for CL4 relative to DOPG (Figure 3). This result is consistent with the decrease in area per lipid, which would lead to chain compaction. For all bonds, DOPG showed an order parameter closer to zero than CL4, which is a consequence of the added motional freedom in DOPG due to the absence of a head group link. The effect is strongest for the PO4 to head glycerol bond, where CL4 shows a strong tendency to align in the membrane plane. In DOPG, this effect is smaller but still shows a preference for the perpendicular orientation. When compared to dipalmitoylphosphatidylcholine (DPPC),36 both DOPG and CL4 have more negative order in the terminal or linking bond, but overall, the order parameters are similar. Phase Polymorphism of CL Derivatives. The main result of our simulations is the phase preference of the CL derivatives at different pHs, as identified by inspecting the aggregate shapes in Figure 4. Some features of the lipid structures in the simulations require further definitions. The intermediate structure between the lamellar and HII phases contain stalks, which are connections between adjacent bilayers. The stalks were previously found to constitute the rate-limiting step of lamellar to hexagonal phase transitions.37 Transition intermediates with stalks from simulations of CL4 at low pH, where the phase transition was completed, were used as starting structures for other simulations, where the head group charge and number of acyl chains was changed before equilibrating the system. Another feature made possible by the special head group in CL is a splayed-chain lipid conformation,53 where acyl chains from different parts of the lipid are inserted into different aggregates, effectively linking them together. Splayed-chain lipids were mainly found when simulations were started from isotropic mixtures, and thus may not represent a thermodynamically stable state. They were, however, persistent for hundreds of nanoseconds. Dilysocardiolipin (CL2). For CL2, two main structures were obtained, the micellar phase at neutral pH and the lamellar phase at low pH (Figure 4, 2A and 2B, respectively). Simulations of CL2 were started from isotropic mixtures obtained from short (1.6 ns) high temperature runs (2000 K) with constant volume. In Figure 4, 2D, a typical cross section of the cylindrical micellar phase is shown. A number of splayed-chain lipids were found between micelles at neutral pH and between bilayers (Figure 4, 2E) at low pH. Many splayed-chain lipids persisted in this conformation the entire duration of the simulations. The transition from the micellar phase to the lamellar phase at low pH was not inhibited by the high water content, indicating that the mechanism proceeds locally as soon as head group charges are neutralized. The process leads to increased micelle size and eventually forms bilayers with water filled pores. As the pores are eliminated, stable bilayers remain. Changing from a rectangular to a triclinic box and restarting the simulations
7198 J. Phys. Chem. B, Vol. 111, No. 25, 2007
Dahlberg
Figure 4. Polymorphism of cardiolipin derivatives. Water and ions have been removed for clarity: (1) CG representation of the lipids. (2A) CL2 in micelles, with the head groups rendered thinner in order to show the hydrocarbon core. The periodic box is shown in black. (2B) Lamellar CL2 with splayed-chain lipids. (2C) Micellar phase with 8000 CL2 lipids. Water and ions are shown in blue. (2D) Cross section of a cylindrical CL2 micelle. (2E) Close-up of 2B with a splayed-chain lipid shaded in blue. (3A) Lamellar CL3 phase. (3B) Lamellar CL3 phase with stalks. (4A) Lamellar CL4 phase. (4B) Inverse hexagonal CL4 phase. (5A) Inverse hexagonal CL5 phase. Note the lighter shaded chain particle at the end of acyl tails. Slightly different scales are used to show the characteristics of each phase.
did not change the final structure, which indicates that the lamellar phase is not an artifact of box geometry. The mechanism of transition is quite different between the lamellar and HII phases, where stalk formation between bilayers constitutes the rate-limiting step.37 Upon decreasing the water content, branched micelles appeared. Similar effects can be seen with changes in lipid geometry, as evidenced by dissipative particle dynamics studies of single chain lipids.54 These phases are consistent with those found by Powell and Marsh. In the absence of salt, the 31P spectrum consisted of a narrow, isotropic peak indicative of the micellar phase. At high salt, evidence for axially symmetric motional averaging was found, supporting the hypothesis of a lamellar phase under these conditions.18 Monolysocardiolipin (CL3). The behavior of CL3 indicates that the lamellar phase is favored for both pHs studied here. The 31P spectra obtained by Powell and Marsh changed very little upon the addition of NaCl and remained in the lamellar phase even at 3 M NaCl.18 In order to assess phase stability, transition intermediates between the lamellar and HII phase were used as starting structures. At neutral pH, the system reverted into the lamellar phase (Figure 4, 3A). At low pH, some stalks disappeared quickly, whereas others remained without elongating for the duration of the simulation (Figure 4, 3B). Splayedchain lipids at the point of contact between the bilayers likely influence the longevity of the stalks. The CL3 simulations were initialized from a lamellar phase, and hence run the risk of being biased toward bilayers rather than the micelles. To avoid this bias, a simulation of a smaller system (96 lipids) at high hydration and neutral pH was started from an isotropic mixture of lipids, ions, and water. A pore-free bilayer assembled rapidly. If a micellar structure would be favorable, it would likely manifest itself when isotropic starting conditions are applied. The resulting bilayer supports the conclusion that the lamellar phase is preferred for CL3 at neutral pH. Cardiolipin (CL4). CL4 bilayer stacks with eight waters per lipid (Figure 4, 4A) were simulated, and no spontaneous stalk formation was seen for 160 ns. When started from a stalk intermediate, constructed from simulations at low pH, the system did not proceed into the HII phase, but stalks remained without elongation for 1.2 µs with a triclinic box. The same result was
obtained for a rectangular geometry in a 160 ns simulation. These results are consistent with CL4 preferring the lamellar phase at neutral pH, which was also found by Powell and Marsh with 31P NMR at low salt content.18 At low pH, bilayer stacks made the transition into the HII phase (Figure 4, 4B) spontaneously. Although the bilayer stack was built with one in-plane dimension significantly shorter than the other, hexagonal structures were formed in either of the two directions parallel to the bilayer plane, with a resulting hexagonal spacing of 6.8 and 6.9 nm, respectively. This result reduces the concern for finite-size effects due to periodic boundary conditions, because it shows that the transition is not dependent on the smallest dimension of the simulation box. Because the coarse-grained model employs five particles in the chains, where ideally 4.5 should be used for an 18-carbon chain, an approximately 10% too large spacing is obtained. The spacing was still within the range 6-7 nm found experimentally by Seddon et al.,14 and the general HII phase preference is supported by 31P NMR.18,55 Although, as noted previously, the ion distribution close to the membrane changes substantially as the charge interaction cutoff was increased from 1.2 to 2.0 nm, we saw the lamellar to HII phase transition also in the CL4 system at low pH with a 1.2 nm cutoff. The transition was even slightly faster (120 and 100 ns for 2.0 and 1.2 nm, respectively), implying that longrange electrostatics play a role in stabilizing the lamellar phase or inhibiting stalk formation and elongation. Acylcardiolipin (CL5). Due to the fifth acyl chain, CL5 does not have a free hydroxyl group on the central glycerol. The effect with respect to the dissociation constants is probably similar to the 2′-deoxy analogue of cardiolipin where both apparent pK values are below 4.56 Hence, CL5 is most likely doubly anionic at neutral pH, which motivates using a model with charges on both PO4 groups. In simulations started from a lamellar structure with eight waters per lipid, stalks appeared quickly, and the formation of a hexagonal phase with many defects ensued. The numerous persistent defects found in the CL5 structure at low hydration are consistent with previous observations of phase transitions in CG lipids.37 This behavior is to be expected when a system is prepared far from equilibrium, because fast and local
Polymorphic Phase Behavior of CL Derivatives
Figure 5. Acyl chain angle (θ) for the lamellar to HII transition in CL5 over the course of 2 µs. The angle distributions are averaged over 200 ns each and begin at 0 µs (red), 0.6 µs (orange), 1.2 µs (green), and 1.8 µs (blue), respectively. Snapshots from the simulation are shown with three periodic boxes (one simulation box outlined at 0 µs in black). A flexible triclinic box was used to allow relaxation during equilibration and phase transition.
transitions can trap the system in metastable conformations, which in this case are locally hexagonal, but with low longrange ordering. In order to slow down the transition into the hexagonal phase, higher water per lipid content (12 waters per lipid) was used with triclinic boundary conditions. In a simulation 2 µs in length, spontaneous formation of a defect-free HII phase was observed (Figure 4, 5A, and Figure 5). At this time scale, there is time for a cooperative transition, slow enough to avoid metastable states with defects. As a measure of the molecular geometry, we measured the angle (θ) between the vectors from the linking glycerol particle to the terminal chain particle at the sn-1 and sn-1′ acyl chains (see Figure 5). Distributions of θ for windows of 0.2 µs length were calculated and are presented in Figure 5. As the transition from the lamellar to the HII phase proceeds, the distribution of angles is widened and shifted toward larger values, reflecting the connection between the molecular geometry and the curvature of the phase. The small persistent peak at 13° corresponds to the minimum in Lennard-Jones energy between the two terminal chain particles. Discussion Experimentally, the phase transitions for cardiolipin derivatives appear at very high salt (3-4 M) or at pH below the first apparent pK value for CL4. At a NaCl concentration of 1 M, Seddon et al.14 found the phase transition when the pH was lowered to 2.8. With corrections for electrostatics,10 the intrinsic second pK value lies between 4 and 6, where the apparent pK values were between 7.5 and 9.5.56 In the simulations presented here, phase transitions were observed between the fully charged model and the half charged model. Mapped to the intrinsic pK values,10 the fully charged model is in the neutral pH range. As half the charge is removed, and phase transitions occur, the pH thus corresponds to values below 4-6. Although higher levels
J. Phys. Chem. B, Vol. 111, No. 25, 2007 7199 of unsaturation in CL were used in the experimental study by Powell and Marsh,18 which probably influences both melting temperature and phase behavior, our CL model still has a high degree of unsaturation. For comparison, tetraoleoyl-CL was previously found to adopt the HII phase at 3.5 M NaCl.55 The identity of the acyl chains varies greatly between different species, and even in different tissues within species, but generally has many unsaturated chains.57 A wide variety of combinations are thus possible for CL with four chains. Although the chiral isomers have been lost in the CG model presented here, the degree of unsaturation is still represented by using oleoyl chains with a kink. The tendency for lipids to adopt the HII phase generally increases with the degree of unsaturation, and to a first approximation, we suggest that only the average number of double bonds per chain is important for phase preferences with the CG model. The fully saturated tetramyristoyl-CL studied recently by Lewis et al.51 predictably prefers the lamellar phase but also has a gel to liquid crystalline phase transition at 40 °C. Although the naturally occurring CLs are generally not completely saturated, the experimental study of Lewis et al. provides a good source for future improvement on the CG model. A number of simplifications are made in the coarse-grained model; thus, it is not straightforward to say if the observed shift relative to the apparent pK values is due to the actual intrinsic pK values being lower, or shortcomings in the model. Dramatic changes in pH and salt concentration not only affect charge screening but also affect the hydration at the interface. Although these effects are only captured with small effective differences in the interaction energy between phosphate and water (5.8 kJ/ mol) and phosphate and sodium ions (5.0 kJ/mol), we can still see relevant phase transitions by simply changing the charge of the phosphate group and correcting with sodium ions to attain charge neutrality. Many details in the hydration and ion binding at the interface are lost in this description, and the charge density is very high compared to the inner mitochondrial membrane, where 9.2 mol % CL is observed. Recently, CG simulations of charged dendrimer insertion into DPPC membranes were reported,58 and the inclusion of longer ranged electrostatic interactions did not appreciably alter the dynamics of the lipids. However, a 1.2 nm cutoff failed to capture the insertion process of dendrimers into the membrane and associated pore formation, indicating that long-range interactions play a serious role when charges are close to the membrane. Although we could reproduce phase transitions for CL4 with a 1.2 nm cutoff, we saw changes in the dynamics (decreased transition time) and support the use of a long cutoff. For the lipids CL3 and CL5, the same water per lipid ratio was used, which changes the relative content of water slightly, due to the removal and addition of acyl chains. The relatively low hydration allows for phase transitions between the lamellar and inverse hexagonal phases in reasonable simulation time. Water content is a strong modifier of phase stability. However, the different hydration levels used in this study were chosen to model phases in which phase separation does not occur, because of the limitations in the length and time scales inherent in the CG MD method. For our CL2 simulations, water content is not as important for the phase transition between the micellar and lamellar phases, but it affects the appearance of the micellar phase. Water per lipid ratios of 18:1 and 37:1 were tried, and resulted in branched and unbranched micelles, respectively. Conclusions The phase preferences derived from the molecular dynamics computer simulations using the coarse-grained models for
7200 J. Phys. Chem. B, Vol. 111, No. 25, 2007 cardiolipin derivatives follow the phase behavior observed experimentally by Powell and Marsh,18 both with respect to (i) the number of acyl chains and (ii) the head group charge (Figure 4). By changing the head group charge directly, as a way of changing the effective pH in the systems, phase transitions between micellar-lamellar and lamellar-inverse hexagonal were observed. Surface charge was found to be an effective modulator of the lipid phase, also reported by, for example, Lewis and McElhaney.59 Interface curvature decreased with decreasing charge and increasing number of oleoyl chains, which is in direct agreement with a change in the head group area and hydrophobic volume in the dimensionless packing parameter.19,20 A further look into the molecular details of cardiolipin charge is warranted, and will be the objective of future investigations. Acknowledgment. The author thanks Profs. Arnold Maliniak and Aatto Laaksonen for comments on the manuscript. References and Notes (1) Gomez, B., Jr.; Robinson, N. C. Anal. Biochem. 1999, 267, 212216. (2) Hovius, R.; Lambrechts, H.; Nicolay, K.; de Kruijff, B. Biochim. Biophys. Acta 1990, 1021, 217-226. (3) Haines, T. H.; Dencher, N. A. FEBS Lett. 2002, 528, 35-39. (4) Hoch, F. L. Biochim. Biophys. Acta 1992, 1113, 71-133. (5) Schlame, M.; Rua, D.; Greenberg, M. L. Lipid Res. 2000, 39, 257288. (6) Gohil, V. M.; Hayes, P.; Matsuyama, S.; Schagger, H.; Schlame, M.; Greenberg, M. L. J. Biol. Chem. 2004, 279, 42612-42618. (7) Pfeiffer, K.; Gohil, V.; Stuart, R. A.; Hunte, C.; Brandt, U.; Greenberg, M. L.; Schagger, H. J. Biol. Chem. 2003, 278, 52873-52880. (8) Zhang, M.; Mileykovskaya, E.; Dowhan, W. J. Biol. Chem. 2005, 280, 29403-29408. (9) McAuley, K. E.; Fyfe, P. K.; Ridge, J. P.; Isaacs, N. W.; Cogdell, R. J.; Jones, M. R. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 14706-14711. (10) Gorbenko, G. P.; Molotkovsky, J. G.; Kinnunen, P. K. Biophys. J. 2006, 90, 4093-4103. (11) Vreken, P.; Valianpour, F.; Nijtmans, L. G.; Grivell, L. A.; Plecko, B.; Wanders, R. J.; Barth, P. G. Biochem. Biophys. Res. Commun. 2000, 279, 378-382. (12) Hauff, K. D.; Hatch, G. M. Prog. Lipid Res. 2006, 45, 91-101. (13) Schlame, M.; Ren, M. FEBS Lett. 2006, 580, 5450-5455. (14) Seddon, J. M.; Kaye, R. D.; Marsh, D. Biochim. Biophys. Acta 1983, 734, 347-352. (15) Allegrini, P. R.; Pluschke, G.; Seelig, J. Biochemistry 1984, 23, 6452-6458. (16) Cullis, P. R.; Verkleij, A. J.; Ververgaert, P. H. Biochim. Biophys. Acta 1978, 513, 11-20. (17) de Kruijff, B.; Cullis, P. R. Biochim. Biophys. Acta 1980, 602, 477490. (18) Powell, G. L.; Marsh, D. Biochemistry 1985, 24, 2902-2908. (19) Cullis, P. R.; de Kruijff, B. Biochim. Biophys. Acta 1979, 559, 399420. (20) Israelachvili, J. N.; Marcelja, S.; Horn, R. G. Q. ReV. Biophys. 1980, 13, 121-200. (21) Egorova, E. M. Colloids Surf., A 1998, 131, 7-18.
Dahlberg (22) Egorova, E. M. Colloids Surf., A 1998, 131, 19-31. (23) Nichols-Smith, S.; Kuhl, T. Colloids Surf., B 2005, 41, 121-127. (24) Siegel, D. P. Biophys. J. 1986, 49, 1171-1183. (25) Siegel, D. P. Biophys. J. 1993, 65, 2124-2140. (26) Lyubartsev, A. P.; Laaksonen, A. Phys. ReV. E 1995, 52, 37303737. (27) Lyubartsev, A. P. Eur. Biophys. J. 2005, 35, 53-61. (28) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. J. Chem. Phys. 2004, 120, 10896-10913. (29) Izvekov, S.; Voth, G. A. J. Phys. Chem. B 2005, 109, 2469-2473. (30) Mu¨ller-Plathe, F. ChemPhysChem 2002, 3, 754-769. (31) Shelley, J. C.; Shelley, M. Y.; Reeder, R. C.; Bandyopadhyay, S.; Klein, M. L. J. Phys. Chem. B 2001, 105, 4464-4470. (32) Shelley, J. C.; Shelley, M. Y.; Reeder, R. C.; Bandyopadhyay, S.; Moore, P. B.; Klein, M. L. J. Phys. Chem. B 2001, 105, 9785-9792. (33) Venturoli, M.; Maddalena Sperotto, M.; Kranenburg, M.; Smit, B. Phys. Rep. 2006, 437, 1-54. (34) Nielsen, S. O.; Lopez, C. F.; Srinivas, G.; Klein, M. L. J. Phys.: Condens. Matter 2004, 16, R481-R512. (35) Smit, B.; Hilbers, P. A. J.; Esselink, K.; Rupert, L. A. M.; van Os, N. M.; Schlijper, A. G. Nature 1990, 348, 624-625. (36) Marrink, S. J.; deVries, A. H.; Mark, A. E. J. Phys. Chem. B 2004, 108, 750-760. (37) Marrink, S. J.; Mark, A. E. Biophys. J. 2004, 87, 3894-3900. (38) Kasson, P. M.; Kelley, N. W.; Singhal, N.; Vrljic, M.; Brunger, A. T.; Pande, V. S. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 11916-11921. (39) Shih, A. Y.; Arkhipov, A.; Freddolino, P. L.; Schulten, K. J. Phys. Chem. B 2006, 110, 3674-3684. (40) Shih, A. Y.; Freddolino, P. L.; Arkhipov, A.; Schulten, K. J. Struct. Biol. 2007, 157, 579-92. (41) Bond, P. J.; Holyoake, J.; Ivetac, A.; Khalid, S.; Sansom, M. S. P. J. Struct. Biol. 2007, 157, 593-605. (42) MacCallum, J. L.; Tieleman, D. P. J. Am. Chem. Soc. 2006, 128, 125-130. (43) Niemela, P. S.; Ollila, S.; Hyvonen, M. T.; Karttunen, M.; Vattulainen, I. PLoS Comput. Biol. 2007, 3, e34. (44) de Vries, A. H.; Yefimov, S.; Mark, A. E.; Marrink, S. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 5392-5396. (45) Knecht, V.; Marrink, S. J. Biophys. J. 2007, 92, 4254-4261. (46) Gurtovenko, A. A.; Miettinen, M.; Karttunen, M.; Vattulainen, I. J. Phys. Chem. B 2005, 109, 21126-21134. (47) Zhao, W.; Rog, T.; Gurtovenko, A. A.; Vattulainen, I.; Karttunen, M. Biophys. J. 2007, 92, 1114-1124. (48) Elmore, D. E. FEBS Lett. 2006, 580, 144-148. (49) Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (50) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306-317. (51) Lewis, R. N.; Zweytick, D.; Pabst, G.; Lohner, K.; McElhaney, R. N. Biophys. J. 2007, 92, 3166-3177. (52) Petrache, H. I.; Tristram-Nagle, S.; Gawrisch, K.; Harries, D.; Parsegian, V. A.; Nagle, J. F. Biophys. J. 2004, 86, 1574-1586. (53) Corkery, R. W. Colloids Surf., B 2002, 26, 3-20. (54) Cooke, I. R.; Deserno, M. Biophys. J. 2006, 91, 487-495. (55) Sankaram, M. B.; Powell, G. L.; Marsh, D. Biochim. Biophys. Acta 1989, 980, 389-392. (56) Kates, M.; Syz, J. Y.; Gosser, D.; Haines, T. H. Lipids 1993, 28, 877-882. (57) Schlame, M.; Ren, M.; Xu, Y.; Greenberg, M. L.; Haller, I. Chem. Phys. Lipids 2005, 138, 38-49. (58) Lee, H.; Larson, R. G. J. Phys. Chem. B 2006, 110, 18204-18211. (59) Lewis, R. N.; McElhaney, R. N. Biophys. J. 2000, 79, 1455-1464.
