2084
Langmuir 2001, 17, 2084-2096
Bending Frustration of Lipid-Water Mesophases Based on Cubic Minimal Surfaces1 U. S. Schwarz†,2 and G. Gompper*,‡ Department of Materials and Interfaces, Weizmann Institute of Science, Rehovot 76100, Israel, and Institut fu¨ r Festko¨ rperforschung, Forschungszentrum Ju¨ lich, 52425 Ju¨ lich, Germany Received October 3, 2000. In Final Form: January 2, 2001 Inverse bicontinuous cubic phases are ubiquitous in lipid-water mixtures and consist of a lipid bilayer forming a cubic minimal surface, thereby dividing space into two cubic networks of water channels. For small hydrocarbon chain lengths, the monolayers can be modeled as parallel surfaces to a minimal midsurface. The bending energy of the cubic phases is determined by the distribution of Gaussian curvature over the minimal midsurfaces which we calculate for seven different structures (G, D, P, I-WP, C(P), S, and F-RD). We show that the free-energy densities of the structures G, D, and P are considerably lower than those of the other investigated structures due to their narrow distribution of Gaussian curvature. The Bonnet transformation between G, D, and P implies that these phases coexist along a triple line, which also includes an excess water phase. Our model includes thermal membrane undulations. Our qualitative predictions remain unchanged when higher order terms in the curvature energy are included. Calculated phase diagrams agree well with the experimental results for 2:1 lauric acid/dilauroyl phosphatidylcholine and water.
I. Introduction When dissolved in aqueous solvent, lipids self-assemble into an amazing variety of different structures as a function of concentration and temperature.4-6 Most prominent is the lamellar phase, which consists of a stack of lipid bilayers separated by layers of water. In fact, the lipid bilayer constitutes the basic building block of biological membranes; in eucaryotic cells, it not only envelopes the cell itself but also separates its nucleus and its organelles from the cytosol and builds up the endoplasmic reticulum, the Golgi apparatus, and the unilamellar vesicles which serve as reaction chambers and transport vehicles. The predominance of the lamellar phase at ambient temperatures derives from the fact that in contrast to surfactants, which usually have large headgroups and form micelles, lipids have rather bulky hydrocarbon chains. A so-called bilayer lipid has a chain region cross section which is of similar size as the area per headgroup; thus the spontaneous curvature is small and the lipid forms a lamellar phase. If temperature is increased, the size of the headgroups’ hydration shell decreases, spontaneous curvature increases, and the lipid monolayers tend to curve toward the water regions. Since the tendency of lipid monolayers to curve is frustrated in the lamellar phase, it is often found to transform into a inverse hexagonal phase which * To whom correspondence may be addressed. † Department of Materials and Interfaces, Weizmann Institute of Science. ‡ Institut fu ¨ r Festko¨rperforschung, Forschungszentrum Ju¨lich. (1) Parts of this publication were already presented in a Letter: Schwarz, U. S.; Gompper, G. Phys. Rev. Lett. 2000, 85, 1472-1475. (2) Present address: Max-Planck-Institut fu¨r Kolloid- und Grenzfla¨chenforschung, 14424 Potsdam, Germany. (3) Abbreviations: IBCP, inverse bicontinuous cubic phase; TPMS, triply periodic minimal surface; IPMS, infinite periodic minimal surface; LA, lauric acid; DLPC, dilauroyl phosphatidylcholine. (4) Fontell, K. Colloid Polym. Sci. 1990, 268, 264-285. (5) Luzzati, V.; Vargas, R.; Mariani, P.; Gulik, A.; Delacroix, H. J. Mol. Biol. 1993, 229, 540-551. (6) Seddon, J. M.; Templer, R. H. Polymorphism of lipid-water systems. In Structure and dynamics of membranes-from cells to vesicles; Lipowsky, R., Sackmann, E., Eds.; Elsevier: Amsterdam, 1995; Vol. 1A.
consists of a two-dimensional array of inverse cylindrical micelles. Spontaneous curvature can also be increased by changing molecular architecture, by adding lipids with bulkier chains, or by increasing salt concentration (since salt screens electrostatic repulsion between charged headgroups).7 Lipids with spontaneous curvature are often called nonbilayer lipids. For large spontaneous curvature, a cubic phase of inverse spherical micelles is sometimes observed.6 The inverse phases have been denoted type II (water-in-oil) phases in contrast to type I (oil-in-water) phases, which dominate surfactant phase behavior.8 Although there are some lipids which also form type I phases, we focus here on type II phases since these are the ones formed by the majority of biologically relevant lipids. One of the intriguing aspects of the polymorphism of lipid-water mixtures is the existence of yet another structural type, which often occurs between the lamellar and the inverse hexagonal phase. Like the lamellar phase, these inverse bicontinuous cubic phases (IBCPs) consist of lipid bilayers. However, now a single lipid bilayer extends throughout the whole of space and divides it into two disconnected but interpenetrating labyrinths of cubic symmetry, which both are filled with water. Thus the structure can also be considered to be built from water channels (that is, inverse cylindrical micelles) as is the inverse hexagonal phase, only that now the channels meet in vertexes with a coordination number related to the topology of the cubic phase. Although the IBCPs unite the structural elements of the lamellar and the inverse hexagonal phases, in contrast to them they are optically isotropic and highly viscous. Like the inverse hexagonal phase, the IBCPs are stabilized by spontaneous curvature. Therefore they are formed by non-bilayer lipids and their stability can be controlled by changing the lipids’ spontaneous curvature. The generic occurrence of IBCPs in lipid-water mixtures and their property of dividing space into two (7) Porte, G. J. Phys.: Condens. Matter 1992, 4, 8649-8670. (8) Luzzati, V.; Tardieu, A.; Gulik-Krzywicki, T.; Rivas, E.; ReissHusson, F. Nature 1968, 220, 485-488.
10.1021/la0013805 CCC: $20.00 © 2001 American Chemical Society Published on Web 03/13/2001
Lipid-Water Mesophases
interwoven aqueous compartments lends itself readily to speculations on their biological relevance, e.g., on their possible occurrence in the endoplasmic reticulum or the Golgi apparatus.5,6 A recent survey of a vast number of transmission electron microscopy (TEM) pictures suggests that IBCPs are indeed ubiquitous in biological systems.9,10 The abundance of non-bilayer lipids in biological systems is a longstanding puzzle, and several studies have explored the possible relationship between lipid polymorphism and biological membrane function (for reviews see refs 1113). In particular, it was shown that certain cells such as Escherichia coli regulate the relative amount of bilayer and non-bilayer lipid in order to control membrane morphology and function of membrane proteins. E. coli mutants lacking the non-bilayer lipid phosphatidylethanolamine were found to be severely impaired in regard to protein transport across the plasma membrane.14 Luzzati and co-workers have pointed out that if any IBCP had a biological function, it should be stable at physiological temperatures and in excess water.5,12 IBCPs are also relevant for biotechnological applications; e.g., they have been used recently as an artificial matrix which enables membrane proteins to crystallize in a three-dimensional array.15 Since IBCPs consist of a lipid bilayer separating two labyrinths of water, is it reasonable to assume that the two monolayers are arranged in a way which is locally symmetrical about their midsurface. In differential geometry, a surface which curves to both sides in the same way is known as a minimal surface, i.e., a surface for which the mean curvature vanishes at each of its points. Therefore the midsurfaces of the IBCP lipid bilayers are often modeled as triply periodic minimal surfaces (TPMS, sometimes also denoted IPMS for infinite periodic minimal surfaces) of cubic symmetry.16,17 Indeed it has been confirmed by a thorough analysis of electron density maps derived from X-ray data that the midsurfaces of these structures can be very close to minimal surfaces.5 Note that the stronger requirement of the monolayers’ global symmetry about the midsurface leads to a certain subclass of cubic TPMS which is called balanced since then the two water labyrinths have to be congruent to each other. Since the seminal work of Schoen,18 it is known that there exists a large variety of different cubic TPMS. Before Schoen’s work, only P, D, and C(P) were known from the 19th century work of Schwarz and his students; then Schoen described G, F-RD, I-WP, O, C-TO, and C(D). More cubic TPMS have been found later by Karcher and Polthier,19 Fischer and Koch,20 and others, but most of them seem to be too complicated to be of physical relevance. In Figures 2 and 3 we depict surfaces corresponding to (9) Landh, T. FEBS Lett. 1995, 369, 13-17. (10) Hyde, S.; Andersson, S.; Larsson, K.; Blum, Z.; Landh, T.; Lidin, S.; Ninham, B. W. The language of shape; Elsevier: Amsterdam, 1997. (11) Epand, R. M. Chem. Phys. Lipids 1996, 81, 101-104. (12) Luzzati, V. Curr. Opin. Struct. Biol. 1997, 7, 661-668. (13) de Kruijff, B. Curr. Opin. Chem. Biol. 1997, 1, 564-569. (14) Rietveld, A. G.; Koorengevel, M. C.; de Kruijff, B. EMBO J. 1995, 14, 5506-5513. (15) Pebay-Peyroula, E.; Rummel, G.; Rosenbusch, J. P.; Landau, E. M. Science 1997, 277, 1676-1681. Rummel, G.; Hardmeyer, A.; Widmer, C.; Chiu, M. L.; Nollert, P.; Locher, K. P.; Pedruzzi, I.; Landau, E. M.; Rosenbusch, J. P. J. Struct. Biol. 1998, 121, 82-91. (16) Dubois-Violette, E., Pansu, B., Eds. International Workshop on Geometry and Interfaces at Aussois, France. J. Phys., Colloq. 1990, 7. (17) Klinowski, J., Mackay, A. L., Eds. Curved surfaces in chemical structure. Philos. Trans. R. Soc. London, Ser. A 1996, 354, 1969-2192. (18) Schoen, A. H. Minimal surfaces, Technical Report D-5541, Washington, DC, 1970. (19) Karcher, H.; Polthier, K. Philos. Trans. R. Soc. London, Ser. A 1996, 354, 2077-2104. (20) Fischer, W.; Koch, E. Philos. Trans. R. Soc. London, Ser. A 1996, 354, 2105-2142.
Langmuir, Vol. 17, No. 7, 2001 2085
the seven different TPMS treated in this work. Until now, only G, D, and P have been established to be realized as IBCP in lipid-water mixtures.4-6 One system for which all three of these phases is stable is 2:1 lauric acid (LA)/ dilauroyl phosphatidylcholine (DLPC) and water (with P coexisting with excess water).21 The only other system for which G, D, and P are known to be stable is didodecylphosphatidylethanolamine and water;22 however, here G exists with hardly no hydration at all, and the regions for D and P are very difficult to access experimentally, thus this system should not be compared with the theoretical results for fully hydrated phases presented here. In all other systems with IBCPs, only one or two of them are stable (for example, for monoolein and water, G and D are stable, with D coexisting with excess water23). In this work we present a model which can explain why only G, D, and P are observed in lipid-water mixtures and why their actual occurrence seems to depend on the specific system under investigation. The phase diagram predicted by our model agrees well with the experimental phase diagram for 2:1 LA/DLPC and water. Exact (Weierstrass) representations are only known for G, D, P, and I-WP.24-27 However, it has been shown in recent studies28,29 that these and other cubic TPMS can be generated as isosurfaces of a density field, which minimizes the free-energy functional of a simple Ginzburg-Landau model for amphiphilic systems.30 Using the Fourier representations of ref 29, it is now possible to investigate also the physical properties of the IBCPs, for whose minimal midsurfaces no exact representations are known. In particular, the question of whether these structures might correspond to stable IBCPs in lipidwater mixtures can now be addressed quantitatively. In this work we study G, D, P, I-WP, S, C(P), and F-RD, of which G, D, P, S, and C(P) are balanced. The most important contribution to the free energy of amphiphilic interfaces is the curvature energy, which dominates all other contributions if the radii of curvature are large compared with molecular length scales.31 For the inverse phases in lipid-water mixtures, it is generally assumed that another important contribution to the free energy is the stretching energy of the hydrocarbon chains.32-34 To attain optimal free energy, the neutral surfaces of the two monolayers should realize at the same time constant mean curvature c0 and constant distance l to the minimal midsurface. Since this is not possible geometrically, the free energy will always be frustrated.35 The relative importance of interface (bending) and bulk (stretching) contributions in self-assembled interfaces is (21) Templer, R. H.; Seddon, J. M.; Warrender, N. A.; Syrykh, A.; Huang, Z.; Winter, R.; Erbes, J. J. Phys. Chem. B 1998, 102, 72517261. Templer, R. H.; Seddon, J. M.; Duesing, P. M.; Winter, R.; Erbes, J. J. Phys. Chem. B 1998, 102, 7262-7271. (22) Seddon, J. M.; Hogan, J. L.; Warrender, N. A.; Pebay-Peyroula, E. Prog. Colloid Polym. Sci. 1990, 81, 189-197. (23) Hyde, S. T.; Andersson, S.; Ericsson, B.; Larsson, K. Z. Kristallogr. 1984, 168, 213-219. (24) Lidin, S.; Hyde, S. T.; Ninham, B. W. J. Phys. (Paris) 1990, 51, 801-813. (25) Fodgen, A.; Hyde, S. T. Acta Crystallogr., Sect. A 1992, 48, 442451 and 575-591. (26) Cvijovic, D.; Klinowski, J. Chem. Phys. Lett. 1994, 226, 93-99. (27) Fodgen, A.; Hyde, S. T. Eur. Phys. J. B 1999, 7, 91-104. (28) Go´zdÄ z´, W.; Hołyst, R. Phys. Rev. E 1996, 54, 1-16. (29) Schwarz, U. S.; Gompper, G. Phys. Rev. E 1999, 59, 5528-5541. (30) Gompper, G.; Schick, M. Phys. Rev. Lett. 1990, 65, 1116-1119. (31) Safran, S. A. Adv. Phys. 1999, 48, 395-448. (32) Anderson, D. M.; Gruner, S. M.; Leibler, S. Proc. Natl. Acad. Sci. U.S.A. 1988, 85, 5364-5368. (33) Dan, N.; Pincus, P.; Safran, S. Langmuir 1993, 9, 2768-2771. (34) Du¨sing, P. M.; Templer, R. H.; Seddon, J. M. Langmuir 1997, 13, 351-359. (35) Charvolin, J.; Sadoc, J. F. J. Phys., Colloq. 1990, 51, 83-96.
2086
Langmuir, Vol. 17, No. 7, 2001
a long-debated issue. For amphiphilic systems, for which the free energy contributions come mainly from the interfaces, the main problem is of a mathematical nature: the overall solution can be investigated only in a model which allows the interfacial positions to adjust themselves freely. In contrast, in diblock copolymer systems bulk contributions are much more important and models based on random polymer coils and interacting monomer densities are more adequate. Self-consistent mean-field theory for diblock-copolymer systems suggests that in this case, interfacial tension and stretching are equally important.36 Recently, this framework has been adapted to treat lipidwater mixtures.37 This approach is complementary to the membrane approach; it is most useful for the investigation of systems, in which the internal structure of the bilayer is not uniform. However, self-consistent field theory for large molecules is not expected to describe amphiphilic systems well at large lattice constants, when the membranes’ internal structure is essentially preserved and the main contributions to the free energy depend on their shape. In particular, only in the membrane approach does it becomes possible to address the role of thermally activated membrane undulations. To study interface-dominated systems, it has been suggested in ref 32 to consider the two classes of surfaces, which can completely relax one of the two relevant contributions. Parallel surfaces to the minimal midsurface realize the optimal chain length l and are frustrated only in their bending energy, while constant-mean-curvature surfaces (CMC-surfaces) realize the given spontaneous curvature but are frustrated in their stretching energy. The mathematical properties of parallel surfaces are wellknown in terms of the properties of the underlying minimal surface; the parallel-surface model has therefore been the main tool for the investigation of IBCPs.21,32,38,39 Up to now, the main difficulty of this approach has been to take into account the variation of Gaussian curvature over the minimal midsurface. In earlier work, the curvature energy was therefore expanded in moments of the distribution of the Gaussian curvature. In ref 32 this was done up to eighth order for D and in ref 39 up to second order for P, D, and G. CMC-surfaces21,32,34 have the advantage that they offer a unified picture for all relevant phasess lamellar, inverse bicontinuous cubic, inverse hexagonal, and inverse micellar cubic phases are modeled by planes, cubic CMC-surfaces, cylinders, and spheres, respectively. However, since no exact representations are known, the stretching frustration for IBCPs has so far been evaluated numerically for D only.32 We will show below that in the case of short chains (or large lattice constants), stretching energy is prohibitively large, and the two monolayers of a IBCP in lipid-water systems can be modeled by parallel surfaces. In the case of long chains (or small lattice constants), we expect stretching and bending to be equally important; this agrees with the results for diblock copolymers, which are in the long chain limit. Our reasoning implies that CMC-surfaces are not a good representation of IBCPs in lipid-water systems for physical reasons. However, for small distances to the minimal midsurface, their mathematical properties are very similar to the ones of parallel surfaces;32,40 thus many results from the parallel surface model will carry (36) Matsen, M. W.; Bates, F. S. Macromolecules 1996, 29, 76417644. (37) Li, X.; Schick, M. J. Chem. Phys. 2000, 112, 6063-6072. (38) Hyde, S. T. J. Phys. Chem. 1989, 93, 1458-1464. (39) Helfrich, W.; Rennschuh, H. J. Phys., Colloq. 1990, 51, 189195. (40) Schwarz, U. S.; Gompper, G. J. Chem. Phys. 2000, 112, 37923802.
Schwarz and Gompper
over to the CMC-model. In general, however, the parallel surface model is more appropriate both for physical and for mathematical reasons. Moreover, there is strong experimental evidence for it, both from a detailed reconstruction of electron densities5 and from analysis of swelling data.41 In contrast to earlier work, we evaluate the curvature energy in the parallel surface model without any expansion and to high numerical accuracy by first calculating the distribution of Gaussian curvature over the minimal midsurface. For G, D, P, and I-WP this can be done from their Weierstrass representations. We also consider the structures S, C(P), and F-RD by calculating the distribution from the representations which we obtained from Ginzburg-Landau theory. Our main result is that the free-energy densities of the structures G, D, and P are considerably lower than those of the other investigated structures due to their narrow distribution of Gaussian curvature over the minimal midsurface. We show that this result persists when thermal membrane undulations and higher order terms in the bending energy are considered and argue that potential additional contributions to the free energy are unlikely to change it. This explains why only G, D, and P have been observed in lipid-water mixtures. In fact, the calculated phase diagram agrees nicely with the experimental one for 2:1 LA/DLPC and water, for which these three phases coexist. We show that due to the existence of a Bonnet transformation between G, D, and P, these phases coexist along a triple line in our model. Simultaneously, P coexists with an excess of water due to a mechanism called emulsification failure. The Bonnet transformation also implies that the free-energy densities of G, D, and P scale as a function of concentration with a universal geometrical quantity, which we term topology index. Since the topology index decreases from G to D to P, the gyroid G is most prominent, followed by smaller regions of stability for D and P at higher water concentrations. Any additional contribution to the free energy, which introduces a new length scale, is expected to change the delicate balance between these three phases; this includes stretching contributions and van der Waals and electrostatic interactions. Since such additional effects are specific to a given experimental system, this result explains qualitatively why usually only one or two of these phases are observed in lipid-water mixtures. The paper proceeds as follows. In section II, we specify the free-energy expression which has to be evaluated for each type of minimal midsurface and discuss some immediate consequences of the model. The calculation of the distribution of Gaussian curvature is explained in section III for P, D, G, and I-WP from their Weierstrass representation. Results for the Gaussian-curvature distributions for S, C(P), and F-RD are obtained in section IV from the representations derived recently from a simple Ginzburg-Landau model. These results are combined in section V to predict phase behavior and to numerically calculate phase diagrams. In section VI, we discuss which physical mechanisms will break the Bonnet symmetry between P, D, and G. Finally, we compare our results with experiments and conclude in section VII. II. The Model The bending energy of one lipid monolayer is described by the Canham-Helfrich Hamiltonian42,43 (41) Engblom, J.; Hyde, S. T. J. Phys. II 1995, 5, 171-190. (42) Canham, P. B. J. Theor. Biol. 1970, 26, 61. (43) Helfrich, W. Z. Naturforsch., C 1973, 28, 693-703.
Lipid-Water Mesophases
Eb )
∫ dAl {2κ(Hl - c0)2 + κjKl}
Langmuir, Vol. 17, No. 7, 2001 2087
(1)
where the integration extends over the monolayer’s neutral surface, which is at a distance l from the minimal midsurface of the bilayer. In the framework of the curvature model, bending and stretching contributions decouple at the neutral surface;31 this allows the areastretching contributions to be disregarded. The neutral surface has been shown experimentally to be located close to the polar-apolar interface.44 In the curvature model (1), the neutral surface is characterized by mean curvature H, Gaussian curvature K, and differential area element dA; the superscript l indicates that the neutral surface is a parallel surface. The model parameters are the spontaneous curvature c0, the bending rigidity κ, and the saddle-splay modulus κj of the monolayers. For lipids, the bending rigidity κ for a monolayer is of the order 10-20 kBT. Usually the saddle-splay modulus κj is assumed to have a small negative value; in fact the curvature model of eq 1 without any additional constraints is well-defined only for -2κ e κj e 0.45 The spontaneous curvature c0 can be considered to result from the mismatch between optimal headgroup area, volume incompressibility of the chains, and optimal chain length. For nonionic lipids, it depends mainly on temperature, which changes the headgroup hydration and therefore optimal headgroup area. Since lipids typically form lamellar and inverse phases, where the monolayers prefer to bend toward their polar sides, we define positive curvature to be toward the water regions. We assume that spontaneous curvature scales with temperature as c0 ∼ (T - Tb), since such a linear relationship seems to hold very well for ternary surfactant systems.46 The stretching energy can be assumed to be harmonic about the average chain length l, so that
Es )
∫ dA ks(L - l)2
(2)
where L is the local chain length. The stretching modulus ks is itself a function of the average chain length l. The relevant scaling law has to be obtained from a microscopic model; here, we model the lipid monolayer as a Gaussian polymer brush grafted to the polar-apolar interface, but more realistic models give similar results.31 The stretching energy per chain reads Es ) (3kBT/2N)(l/aK)2, where aK is the (microscopic) Kuhn length and N the number of Kuhn segments. Enforcing the volume constraint ld2 ) NaK3 (where 1/d2 is the grafting density) and transforming from energy per chain to energy per area results in a potential of the form eq 2 with ks ∼ kBT(aK/ld4). Thus the stretching modulus ks scales inversely with average chain length lsstretching becomes very difficult for short chains.31 Although the bending energy eq 1 is well known not to depend on microscopic details, for our purpose it is important to note that the bending rigidity κ also depends on the average chain length l. Starting from the stretching energy of eq 2 and taking into account how the chain volume varies with curvature results in the bending energy of eq 1 with κ ∼ ksl4 ∼ kBT(aKl3/d4).31 The result that the bending rigidity κ scales with the third power of film thickness l is well known also from the elasticity of thin solid sheets. (44) Templer, R. H. Langmuir 1995, 11, 334-340. (45) Helfrich, W. Amphiphilic mesophases made of defects. In Physics of Defects; Balian, R., Kleman, M., Poirier, J.-P., Eds.; Les Houches Summer School 1980; North-Holland: Amsterdam, 1981. (46) Strey, R. Colloid Polym. Sci. 1994, 272, 1005-1019.
We now can estimate the relative importance of bending and stretching contributions. It is known that for both CMC-surfaces and parallel surfaces, the normalized standard deviations of L and H, respectively, are almost independent of hydrocarbon volume fraction,32 so that 〈(L - l)2〉 ∼ l2 and 〈(H - 〈H〉)2〉 ∼ 〈H〉2 ∼ a-2, where a is the lattice constant. Therefore Es/Eb ∼ ksl2/κa-2 ∼ (a/l)2. Thus, for small chain length l (large lattice constant a), the stretching term is much larger and the chains can be assumed to be of constant length. This corresponds to the well-known result that the curvature energy is the relevant contribution to the free energy as long as the curvature radii (that is the lattice constants) are much larger than any molecular length.31 Only for large chain length l (small lattice constant a), the bending energy becomes comparable and both terms have to be considered. Experimentally, this corresponds to the fact that for small water content (i.e., chain length l of comparable size as lattice constant a), the bicontinuous cubic phases become unstable with respect to the hexagonal phase, whose geometrical properties require larger chain stretching. Note that the relative importance of stretching and bending is asymmetric: since the chain length l cannot be larger than the lattice constant a, there is no regime in which bending can be neglected. Thus CMC-surface will never be a good approximation for monolayers in lipid-water systems. However, the stretching frustration can be relaxed by swelling the bilayer with oil; therefore CMCsurfaces are good approximations for amphiphilic monolayers in ternary systems with oil, water, and amphiphile.40 Similar considerations of the relative importance of stretching and bending can be found also in the context of large membrane inclusions (such as transmembrane proteins) where chains have to stretch in order to decrease the hydrophobic mismatch at the inclusion boundary.33 For the rest of this paper, we proceed within the framework of the parallel surface model. We assume that the interfaces between the tail and headgroup regions of the two lipid monolayers in a ICBP are located at distances (l away from a given cubic TPMS; typically l has values between 1 and 2 nm. Note that for a TPMS the mean curvature satisfies H ) (c1 + c2)/2 ) 0 everywhere (where c1 and c2 are the two principal curvatures), while the Gaussian curvature K ) c1c2 ) -c12 varies as a function of the position and is only restricted to satisfy K e 0. The differential area element dAl, the mean curvature Hl, and the Gaussian curvature Kl on the parallel surfaces then follow as functions of l and the quantities dA and K on the minimal midsurface,38,47 with
dAl ) dA(1 + Kl2) Hl )
-Kl 1 + Kl2
Kl )
K 1 + Kl2
(3)
Note that these formulas are special cases of Steiner’s theorem of integral geometry; they are therefore no approximations for small l, but exact. Since minimal surfaces have K e 0, positive mean curvature is defined here to correspond to positive l. The effective bending energy of the lipid bilayer is only a function of the Gaussian curvature K of the minimal midsurface since its mean curvature H vanishes. In an (47) Spivak, M. A comprehensive introduction to differential geometry I-IV; Publish or Perish: Houston, TX, 1979.
2088
Langmuir, Vol. 17, No. 7, 2001
Schwarz and Gompper
Table 1. Euler Characteristic χ, Scaled Surface Area A* (Both in the Conventional Unit Cell), and Topology Index Γ ) (A*3/2π|χ|)1/2 for All TPMS Considered, As Obtained from Local Minima of the Free-Energy Functional of a Simple Ginzburg-Landau Theorya
χ A* Γ
S
G
D
I-WP
P
C(P)
F-RD
-40 5.41457 0.794735
-8 3.09140 0.76665
-16 3.83755 0.74978
-12 3.46367 0.74238
-4 2.34516 0.71637
-16 3.74820 0.655994
-40 4.77522 0.654174
a Exact values are known for G, D, I-WP, P, and C(P) and given in Table 2. The structures are ordered according to decreasing Γ; the value for S is probably too large, but difficult to improve numerically.
expansion for small chain length l, it can be expressed as an infinite series in powers of K. Using eq 1 and eq 3 yields up to order l2
Eb )
∫ dA {4c02l2κ + (2κj + 8c0lκ + 4c02l2κ)K + 4kl2K2} (4)
Therefore the effective saddle-splay modulus for the lipid bilayer, κjbi ) 2κj + 8c0lκ + 4c02l2κ, is not just the sum of the monolayer’s moduli but is corrected to higher positive values due to the presence of the spontaneous curvature c0.7 We conclude that as long as c0l J -κj/4κ, the IBCPs are favored over the lamellar phase, since the preferred curvature of the monolayers translates into a topological advantage of saddle-type bilayer structures. The third term in eq 4 favors saddle-type structures in any case and also removes the degeneracy which minimal surfaces might experience under deformations which preserve the minimal surface property. This has been noted in a seminal work by Bruinsma;48 in his notation, we have κj ) 4κl2 for the elastic modulus of the K2 term. The volume fraction occupied by the hydrocarbon can be calculated as
v)
1 a3
l χ ∫-ll dl′ ∫ dAl′ ) 2A*(al ) + 4π 3 (a)
3
(5)
where we have used the Gauss-Bonnet theorem ∫ dAK ) 2πχ. A* and χ denote surface area and Euler characteristic of the given cubic TPMS in the conventional unit cell with unit lattice constant. For the following, it is useful to introduce the topology index Γ ) (A*3/2π|χ|)1/2. This quantity is independent of scaling and choice of unit cell and characterizes a given TPMS topology in a universal way; the higher its value, the smaller is the porosity and the larger the specific surface area of a TPMS. It is the only quantity which characterizes a two-dimensional minimal surface in three-dimensional space independently of lattice constant and choice of unit cell, and its relevance for the description of bicontinuous cubic phases in amphiphilic systems has been discussed before.27,29,38,40,49 The topology index of the various TPMS is of order 1 and is given in Table 1. The gyroid G is expected to have the highest value since it divides space into labyrinths with 3-fold coordinated vertexes. S has a similar value for Γ, since it features a mixture of 3- and 4-fold coordinated vertexes. Then comes the double diamond structure with 4-fold coordinated vertexes, while all other structures have higher values. In Table 1 we collect the values of A*, χ, and Γ for the cubic TPMS investigated in this work as obtained from our numerical representations in the framework of a simple Ginzburg-Landau theory (see below). Exact values are known for G, D, I-WP, P, and C(P) and given in Table 1. Comparing numerical and exact (48) Bruinsma, R. J. Phys. II 1992, 2, 425-451. (49) Stro¨m, P.; Anderson, D. M. Langmuir 1992, 8, 691-709.
Figure 1. Lattice constant a in units of the chain length l as a function of hydrocarbon volume v. We use a/l ) v/2A*, which is an excellent approximation up to v ≈ 0.8, beyond which the inverse bicontinuous cubic phases as modeled here cannot exist anymore due to self-intersection. From left to right, the curves correspond to P, G, I-WP, C(P), D, F-RD, and S, respectively (the curves for I-WP and C(P) nearly collapse and therefore appear as an apparent bold line).
values shows that for most structures, our numerical results are quite accurate. Due to numerical limitations for the more complicated structures, the values for S and F-RD are less precise. In particular, we expect S to have a smaller value for Γ than for G. Equation 5 can be inverted numerically to give the dimensionless function a/l, the lattice constant a in units of the chain length l, as a function of hydrocarbon volume v. For small v we find
(
)
1 a 2A* ) 1+ v2 + O(v4) 2 l v 12Γ
(6)
In Figure 1 we plot the first-order approximation a/l ) 2 A*/v of eq 6 for all structures considered; it can be hardly distinguished from the full curve over the full range of v. By visual inspection of our numerical representations (see below), we find that for all structures considered, parallel surfaces do not self-intersect as long as v j 0.8; IBCPs can accommodate large amounts of hydrocarbon since their geometry is so close to that of the lamellar phase. A typical value for l is 15 Å, so that for v ) 0.5 we have aS ) 325 Å, aF-RD ) 285 Å, aD ) 230 Å, aCP ) 211 Å, aI-WP ) 208 Å, aG ) 185 Å, and aP ) 141 Å. Note that at a given v, the ratio of the lattice constants of two coexisting IBCPs is simply the ratio between their scaled surface areas A*. In Figures 2 and 3 we depict for each of the structures considered one of the two monolayers for the hydrocarbon volume fraction v ) 0.5. Then the volume fraction of each of the two labyrinths is (1 - v)/2 ) 0.25. We are now in a situation to rewrite the curvature energy from eq 1 for our purpose. We consider a typical experimental situation, in which temperature is controlled at constant volume. As explained above, changing temperature T amounts to changing spontaneous curvature c0. Under the assumption of incompressibility of both
Lipid-Water Mesophases
Langmuir, Vol. 17, No. 7, 2001 2089
coexistence, it is not clear how it will distribute into the two phases. The controlled quantity can be considered to be the chemical potential coupled to v (that is the derivative of free energysor free energy volume densitysfor hydrocarbon volume fractionsor interfacial area), and therefore phase transitions have to be calculated by using the Maxwell construction for the free energy density as a function of v (which amounts to requiring the same chemical potential for both phases). We want to remark parenthetically that the free energy per unit area cannot be used to calculate phase coexistences, since an exchange of area, controlled by the chemical potential of the lipid, has to be allowed between coexisting phases. We write the free energy volume density dimensionless by using a factor l/(2κc02). For the lamellar phase the hydrocarbon volume fraction v ) 2 l/a, thus its dimensionless free-energy density follows from eq 1 to be
fL )
Figure 2. One of the two monolayers of P, D, G and I-WP for hydrocarbon volume fraction v ) 0.5. These structures have space groups Pm3 h m, Fd3 h m, I4132, and Im3 h m, respectively. In the full structures, a second monolayer exists parallel to the one shown. Each monolayer defines a space-percolating network filled with water. For the balanced cases P, D, and G, the second network is congruent to the first; this additional symmetry changes the structures’ space groups to Im3 h m, Pn3 h m, and Ia3 h d, respectively. For I-WP, the two networks are different and only the monolayer corresponding to the I network is shown.
1 l 4κc02A ) v 2 Aa 2κc0
(7)
The larger v, the smaller a and the more frustrated curvature energy per volume accumulates. For IBCPs, we again use eq 1 and eq 3 in order to write the dimensionless free-energy density as a function of hydrocarbon volume fraction v and the properties of the bilayers minimal midsurface
( )∫
l fb ) 2 a
(
(
)
(al ) + 1 ( ) ) 1 + K* l c l ( (a) ) 4πχ l r ( ) (8) c l a
l dA* 1 + K* a
2
2
2
K*
2
0
3
2 2
0
where again we have used the Gauss-Bonnet theorem and defined r ) -κj/2κ. As mentioned above, the curvature model eq 1 without any additional constraints is welldefined for 0 e r e 1 and bicontinuous phases are favored for c0l J r/2. A*, K*, and χ refer to a conventional unit cell. Using eq 6 to lowest order, a/l ) 2A*/v, and defining Ξ(K*) ) K*A*/8πχ, we rewrite eq 8 as
{
fb ) v
(
2 -1
v 1 - Ξ(K*)( ) ) ∫ dA* A* ( Γ
1-
( )) v Γ
Figure 3. One of the two monolayers of S, F-RD, and C(P) for a hydrocarbon volume fraction v ) 0.5. These structures have space groups I4 h 3d, Fm3 h m, and Pm3 h m, respectively. In the full structures, a second monolayer exists parallel to the one shown. For the balanced cases S and C(P), the second network is congruent to the first; this additional symmetry changes the structures’ space groups to Im3 h m and Ia3 h d, respectively. For F-RD, the two networks are different and only the monolayer corresponding to the F network is shown.
components, there is only one degree of freedom for composition, which we take to be the hydrocarbon volume v (which is equivalent to the amount of interfacial area). Dividing the free energy by the constant volume, we arrive at the free energy per unit volume as a function of v. However, v is not a controlled quantity, since at two-phase
1 + c0 Ξ(K*) × c0 2 2
+
( ) } (9)
r v 4c02 Γ
2
In eq 9 and below, spontaneous curvature c0 is measured in units of 1/l. Since for the minimal midsurface both K* and χ are negative, Ξ(K*) is a positive quantity which varies over the surface. Note that ∫(dA*/A*)Ξ(K*) ) 1/4 due to the Gauss-Bonnet theorem. The free-energy density fL of the lamellar phase, eq 7, is a special case of fb, with Ξ(K*) ) χ ) 0. In this work, we will treat the free-energy density eq 9 without any further approximations. Previously, similar expressions have been expanded in small v (which is equivalent to small l),32,39 similarly as discussed in eq 4. In fact, each term (Kl2)n in an expansion for small l corresponds to a term v2n+1 in an expansion for small hydrocarbon volume v. Since Kl2 ) -Ξ(K*)(v/Γ)2, higher-order terms in eq 1 would result in a free-energy density which has a similar structure as eq 8. As we will see below, this means that our results will remain valid in the case that higher-order terms are included in the curvature energy of the monolayers.
2090
Langmuir, Vol. 17, No. 7, 2001
Schwarz and Gompper
It is well known that thermal fluctuations can contribute significantly to the free energy of amphiphilic systems. Since we consider the regime where the stretching energy dominates, undulations are the favored fluctuation modes. For the lamellar phase, they give rise to steric repulsion,50 which in our case can be written as
fsteric )
c∞
( )
32c02
kBT 2 v3 κ (1 - v)2
(10)
where c∞ ) 0.106 from Monte Carlo simulations51,52 and field-theoretical calculations.53 For an IBCP, steric repulsion is hardly relevant since the lateral correlation length ξ| is of the order of the lattice constant a. Therefore, the perpendicular correlation length ξ⊥ ∼ (kBT/κ)1/2ξ| , ξ| ∼ a, and few membrane collisions should occur. However, fluctuations of the lipid bilayer lead to a renormalization of κ and κj54 for all phases. Since κ multiplies the average mean curvature squared, which vanishes for the midplane, only the renormalization of κj at length scale l has to be taken into account. We identify the typical length scale of a cubic structure with 〈K〉-1/2. This implies l/l ) 2Γ/v, so that r in eq 9 gets renormalized, with
[
rR ) r -
( )]
2Γ 5 kBT ln 12π κ v
(11)
Since the renormalization increases κj, this effect increases the topological advantage of the saddle-shaped structures as does the spontaneous curvature. Note that for the lamellar phase, these effects might favor the formation of wormholes.55 Here, we neglect this aspect, as well as the acoustic modes of the bicontinuous structures at large wavelengths, which have been discussed by Bruinsma.48 Note also that the thermal fluctuations give a contribution to the free energy, which is linear in the membrane area. We have omitted this term here, since it can be absorbed into the chemical potential of the lipid. Having defined the model, we now turn to its phase behavior. We first note that the dimensionless mean curvature averaged over the parallel surface follows from eq 3 as
∫ dAl Hll (v/Γ)2 ) 〈H 〉ll ) ∫ dAl 4 - (v/Γ)2 l
(12)
where again l/a ) v/2A* and the Gauss-Bonnet theorem has been used. It is plotted in Figure 6a. Note that all curves fall on a universal curve when the hydrocarbon volume fraction v is scaled with the topology index Γ. It reaches the value c0 at
( )
4c0 v) 1 + c0
1/2
Γ
(13)
Thus, c0 is realized in a sequence given by increasing topology index; in particular, the gyroid structure G reaches c0 at the highest value of v. We conclude that if (50) Helfrich, W. Z. Naturforsch., A 1978, 33, 305. (51) Gompper, G.; Kroll, D. Europhys. Lett. 1989, 9, 59-64. (52) Janke, W.; Kleinert, H.; Meinhard, M. Phys. Lett. B 1989, 217, 525-529. (53) Bachmann, M.; Kleinert, H.; Pelster, A. Phys. Lett. A 1999, 261, 127-133. (54) Nelson, D., Piran, T., Weinberg, S., Eds. Statistical mechanics of membranes and surfaces. In Jerusalem winter school for theoretical physics; World Scientific: Singapore, 1989; Vol. 5. (55) Gompper, G.; Goos, J. J. Phys. II 1995, 5, 621-634.
all considered structures were stable, they would appear in the sequence G-S-D-I-WP-P, etc., with increasing water concentration. The stability of the different phases is determined by the free-energy density given in eq 8. The second (topological) term is easy to understand. The topological properties of the parallel and midsurfaces are the same: for negative saddle-splay modulus κj (r > 0), the different IBCPs are favored according to their values of the topology index Γ. Here the G structure performs best since it has the lowest porosity. The first term in eq 8 is more complicated, since the mean-curvature properties of the parallel surfaces translate into the Gaussian-curvature properties of the minimal midsurface in a quite complicated way. As discussed above, the spontaneous curvature of the monolayers implies a topological advantage of the bicontinuous phases. For a more detailed analysis, which includes all terms of the expansion, we note that the first term in eq 8 measures the standard deviation of the parallel surface’s mean curvature from the spontaneous curvature. Since it follows from eq 3 that the distribution of mean curvature over the parallel surface, HldAl, is proportional to the distribution of Gaussian curvature over the minimal surface, KdA, we expect that those structures will be more favorable which have small standard deviations for their Gaussian curvature distributions.39 To quantify this concept, we now turn to the distributions of Gaussian curvature which later will allow us to evaluate eq 8 without any further approximations. III. Curvature Properties from Weierstrass Representations Weierstrass representations are known for P, D, and G25,27 and I-WP.24,26 For each of these TPMS, a fundamental domain can be identified, so that the rest of the surface follows by replicating it with the appropriate space group symmetries (Im3 h m, Pn3 h m, Ia3 h d, and Im3 h m, respectively). The Weierstrass representation is a conformal mapping of certain complicated regions within the complex plane onto the fundamental domain
(x1,x2,x3) ) Re
∫0u+iv dz R(z)(1 - z2, i(1 + z2), 2z)
(14)
where (u, v) are the internal (and conformal) coordinates of the minimal surface. The geometrical properties of such a surface follow as
dA(z) ) |R(z)|2(1 + |z|2)2 du dv H(z) ) 0 K(z) )
-4 |R(z)|2(1 + |z|2)4
(15)
with z ) u + iv. Obviously the (isolated) poles of R(z) correspond to the flat points (K ) 0) of the minimal surface. Only few choices of R(z) yield embedded minimal surfaces. The ones for D and P have been known since the 19th century from the work of Schwarz: for D it is R(z) ) (z8 - 14z4 + 1)-1/2. P corresponds to the same region and follows simply by the Bonnet transformation R(z) f eiθR(z) with θ ) 90° Equation 15 implies that P and D have the same metric and the same distribution of Gaussian curvature. However, since they map differently into embedding space, they have different space groups and lattice constants. The gyroid G was discovered in 1970 by Schoen18 as another Bonnet transformation of D, with θ ) 38.015°. The Weierstrass representation for I-WP was
Lipid-Water Mesophases
Langmuir, Vol. 17, No. 7, 2001 2091
Table 2. Scaled Surface Area A*, Euler Characteristic χ in the Conventional Unit Cell, and Topology Index Γ ) (A*3/ 2π|χ|)1/2 for Those TPMS, for Which Exact Results Are Availablea G D I-WP P C(P)
χ
A*
Γ
N
a
-8 -16 -12 -4 -16
3(1 + k22)/2k2 ) 3.091444 3/k2 ) 3.837785 2(31/2) ) 3.464102 3k2 ) 2.345103 3/k3 ) 3.510478
0.766668 0.749844 0.742515 0.716346 0.655993
24 48 48 12
8k1/3(1 + k22)1/2 ) 2.656243 8k1/3 ) 3.37150 3(31/2)k43/4π ) 7.949874 4k1/3k2 ) 2.156516
a For surfaces with known Weierstrass representations, we also give N, the number of fundamental domains needed to build up the surface in the conventional unit cell, and a, the resulting lattice constant. Here k1 ) F(31/2/2, 81/2/3) where F is the incomplete elliptic integral of the first kind; k2 ) K(1/2)/K(31/2/2) where K(k) ) F(1,k) is the complete elliptic integral of the first kind; k3 ) K(1/31/2)/K(2/3)1/2; and k4 ) Γ(1/3). The surface areas and Euler characteristics of P, D, and G are related to each other due to their Bonnet transformations, as described in the text. The Bonnet angle θ ) 38.015° for G follows from tan θ ) k2. Note that often a smaller unit cell is chosen for D; then one has χ ) -2 and A* ) 1.9188925.
Figure 4. Fundamental domains for Weierstrass representations: (a) G, D, and P; (b) I-WP. Filled circles mark the poles of the generating functions R(z) (which correspond to the flat points of the surfaces) and the hatched regions correspond to the fundamental domains of the TPMS.
found only recently.24,26 If one of its poles is chosen to be at infinity, one has R(z) ) (z(z4 + 1))-2/3. In Figure 4 we show the fundamental domains in the complex plane of the different Weierstrass representations. In Table 2 we collect the values for scaled surface area A* and lattice constant a as they follow from the Weierstrass representations. We also give the number N of replications of the fundamental domain needed to build up the surface in one conventional unit cell and the values for the Euler characteristic χ which follows from the topology generated in this process. Note that for D one can choose another conventional unit cell which is contained n ) 8 times in the one chosen here. This has been done in refs 21 and 32, then χ, A, and K scale as n-1, n-1/3, and n-2/3, respectively. In particular, one then has χ ) -2 and A ) 1.918 892 5. Our choice is motivated by the fact that it is the relevant one also for ternary systems. We define the distribution function of Gaussian curvature
f(K) )
∫ dA(u, v) δ(K - K(u,v))
Figure 5. Distributions of Gaussian curvature for (a) I-WP, D, G, and P as obtained from the Weierstrass representations and (b) S, C(P), and F-RD as obtained from our numerical representations as isosurfaces to Ginzburg-Landau local minima. Since P, D, and G are related by a Bonnet transformation, they differ only by simple rescaling. These data are also shown in ref 29.
By this procedure, we have
fi )
+∆K/2 dK′ δ(K′ - K(u,v)) ) ∆K f(Ki) ∫ dA(u,v) ∫KK-∆K/2 i
i
(17)
(16)
over some surface parametrized by internal coordinates (u, v), where in our case dA(u,v) and K(u,v) are determined by eq 15. For the generating functions R(u,v) detailed above, it is not possible to calculate f analytically. Therefore, we evaluate it numerically by considering a histogram {fi} for a set of discrete Gaussian curvature values {Ki} spaced equidistantly with ∆K. We cover the (u,v) region corresponding to the fundamental domain with a square grid of M × M points. For each point, we calculate dA and K from eq 15 and add dA to the fi which corresponds to the Ki with K ∈[Ki - ∆K/2, Ki + ∆K/2]. To obtain results for one conventional unit cell, the values {fi} are then multiplied by the number N of replications and the values {Ki} are scaled with a2 (both N and a are given in Table 2). In practice we use M ) 30 000 and ∆K ) 0.025; the latter value amounts to about 150 bins for the histogram.
The moments SN of the distribution f(K) follow as
SN )
∫ dK f(K)KN ) ∫ dA KN ) ∑fiKiN
(18)
i
In particular, S0 ) ∑ifi ) A* and S1 ) ∑ifiKi ) 2πχ. The procedure described here was carried out to obtain {Ki,fi} for D and I-WP. Since P and G are related to D by a Bonnet transformation, their distributions can be obtained simply by appropriately rescaling. For example, for the transformation D f G, the Ki values have to be rescaled with aG2/aD2 and the fi with AG*/AD*, as given in Table 2. In Figure 5a, we plot the four different distributions obtained. P, D, and G feature one peak around a (relatively large) K value, which means they are quite uniformly curved, although they have flat regions, too
2092
Langmuir, Vol. 17, No. 7, 2001
Schwarz and Gompper
Table 3. First Moments SN and Variance ∆ ) S0S2/S12 - 1 of the Distributions of Gaussian Curvaturea S0 S1 S2 S3 S4 S5 S6 S7 S8 ∆
D
I-WP
S
F-RD
C(P)
3.83763 -100.531 3209.46 -110842 4.00533 × 106 -1.49309 × 108 5.69880 × 109 -2.21675 × 1011 8.75987 × 1012 0.218702
3.46250 -75.4075 2434.91 -86861.5 3.24503 × 106 -1.24566 × 108 4.86942 × 109 -1.92889 × 1011 7.719305 × 1012 0.482666
5.41457 -251.403 18514.0 -1.7220 × 106 1.82465 × 108 -2.08738 × 1010 2.50549 × 1014 -3.10619 × 1014 3.941396 × 1016 0.586079
4.77522 -251.316 21821.3 -2.29671 × 106 2.66794 × 108 -3.27706 × 1010 4.16686 × 1012 -5.42345 × 1014 7.179904 × 1016 0.649801
3.74820 -100.992 5012.37 -286143 1.81085 × 107 -1.21507 × 109 8.45482 × 1010 -6.02391 × 1010 4.363058 × 1014 0.842022
a In particular, S ) A* and S ) 2πχ with A* and χ given in Table 2. If the moments for D are scaled appropriately for a change of unit 0 1 cell, the values given in ref 32 are recovered. The moments for G and P are not shown since they follow from those for D due to the Bonnet transformation. The values for ∆ are the same for G, D, and P.
(which correspond to the poles of R(z)). I-WP is bimodal; that is, it has prominently both flat and strongly curved regions. The strong curvature of D is somehow artificial due to our choice of unit cell. If one corrects for this by a factor 1/4 for the K values, it becomes clear that I-WP has a much broader distribution than P, D, and G. This can also be inferred from Table 3, where we list the first moments SN of the distributions for D and I-WP. The ones for D are identical to the ones given in ref 32 after appropriate rescaling SN with 8-(2N+1)/3 due to the different choice in unit cell (the numerical deviation is less than 0.02% even for N ) 8). The moments for P and G can be calculated from the ones for D due to the Bonnet transformation between them; e.g. SNG ) SND(NG/ND)(aG/ aD)2(N-1) with the values for N and a given in Table 2. In particular, AG* ) AD*(NGaD2)/(NDaG2), χG ) χD(NG/ND), and ΓG ) ΓD(NGaD3)/(NDaG3)sas can be verified from Table 2. It follows immediately that the quantity Ξ(K*) ) K*A*/ 8πχ defined for eq 8 is invariant under a Bonnet transformation. To quantify the variance of the different distributions, we calculate the standard deviation, ∆ ) 〈(K - 〈K〉)2〉/〈K〉2 ) S0S2/S12 - 1 (here 〈. . .〉 means area average). This quantity is independent of scaling, choice of unit cell, and Bonnet transformation and is given in the last row of Table 3. As expected, its value for I-WP is higher than that for P/D/G. IV. Curvature Properties from Ginzburg-Landau To investigate the curvature properties of S, C(P), and F-RD, we use the representations as isosurfaces to fields Φ(r) which we obtained recently as local minima of a simple Ginzburg-Landau theory for ternary amphiphilic systems.29 The same model can also be employed for binary systems if the scalar order parameter field Φ(r) is interpreted to be the local concentration difference between water in the first labyrinth and water in the second labyrinth. Then the midsurface of the lipid bilayer is identified with the surface Φ(r) ) 0. The free-energy functional is given by
F[Φ] )
∫ dr {(∆Φ)2 + g(Φ)(∇Φ)2 + f(Φ)}
(19)
gradients. For balanced TPMS, an additional black and white symmetry has to be implemented. The resulting fields Φ(r) correspond to the various local minima. We showed in ref 29 that by tuning (g0, g2, f0) appropriately, the Φ(r) ) 0 isosurfaces can be made to become very close to minimal surfaces, and generated representations for nine different TPMS: G, D, P, S, C(P), C(D), C(Y), I-WP, and F-RD. K distribution can be obtained from these representations as follows. By triangulating the isosurfaces with the marching cube algorithm, we numerically obtain a metric on the surface. From the Fourier representations of the fields Φ(r), it is possible to calculate exactly the K values on the vertexes of the triangulation. Using a histogram technique like that in the Weierstrass case, it is thus possible to obtain {Ki, fi} even for structures such as S, C(P), and F-RD, which might be of physical relevance but for which no exact representations are known. The distributions obtained in this way for g0 ) -3, g2 ) 7.01, and f0 ) 0 are plotted in Figure 5b. All three of them are roughly trimodal and feature much stronger curved parts than the ones discussed above. This reflects their complicated shape, which makes them more difficult to access numerically than the ones which have been generated by Weierstrass representations; it also explains the difficulty to obtain smooth data for the K distributions. In Table 3, we give their first moments SN and the measure ∆ of their variance. All of them have broader distributions than I-WP, in the sequence S-F-RD-C(P). Therefore we can conclude already from the distributions that the IBCPs should become less favorable in the sequence P/D/G-IWP-S-F-RD-C(P). V. Phase Behavior Given the distribution of Gaussian curvature as histogram {Ki, fi}, it is straightforward to evaluate the first term of the curvature energy, eq 9, of the bicontinuous phases. Since it is an area average over a complicated function of K*, we have to replace ∫ dA* and K* by ∑fi and Ki, respectively. In Figure 6b, we plot the quantity 〈(Hl - c0)2〉ll2 as a function of hydrocarbon volume fraction v for c0 ) 1/6. In contrast to 〈Hl〉ll, compare eq 12, no simple formula exists, and it has to be calculated numerically as
A reasonable choice for f and g has been found to be
f(Φ) ) (Φ + 1)2(Φ - 1)(Φ2 + f0) g(Φ) ) g0 + g2Φ2
(20)
For given model parameters (g0, g2, f0), the free-energy functional is minimized in Fourier space by implementing the correct space group symmetry and minimizing for Fourier amplitudes and lattice constant with conjugate
〈(H - c0) 〉ll ) l
2
2
( )
() ( ( )) () ∑ ( ( )) v
∑ifi 1 - Ξi Γ if i
2
Ξi
v
Γ
1 - Ξi
1 - Ξi
v
Γ
2
2
v
2
+ c0
Γ
2
(21)
Lipid-Water Mesophases
Langmuir, Vol. 17, No. 7, 2001 2093
Figure 7. Free energy densities as a function of hydrocarbon volume fraction v for c0l ) 1/6, r ) 0, and κ/kBT ) 10. The solid line on the right is fL, the other solid lines are the different fb. The IBCPs can be identified from their sequence at the top of the figure, where G, D, P, S, I-WP, C(P), and F-RD appear from right to left. With decreasing v, the phases LR, G, D, and P are stable. The lower dashed line is the Maxwell construction with an excess water phase (emulsification failure). G, D, and P lie on a triple line (compare inset) since they are related by a Bonnet transformation. The right dashed line is the Maxwell construction between G and LR. Reproduced from ref 1. Copyright 2000 American Institute of Physics. Figure 6. (a) 〈Hl〉ll ) (v/Γ)2/(4 - (v/Γ)2) as a function of hydrocarbon volume fraction v. The straight line is c0l ) 1/6. From right to left, the IBCPs appear in the sequence S, G, D, I-WP, P, C(P), and F-RD. (b) 〈(Hl - c0)2〉ll2 as a function of hydrocarbon volume fraction v with c0l ) 1/6. The IBCPs can be identified from their sequence at the top of the figure, where G, D, P, S, I-WP, C(P), and F-RD appear from right to left. The straight line is (c0l)2 ) (1/6)2. Since P, D, and G are related by a Bonnet transformation, their curves differ from each other only by a simple rescaling of v with Γ. Therefore they reach the same minimal values, but at different v according to their Γ values.
Since Ξi is invariant under a Bonnet transformation, 〈(Hl - c0)2〉ll2 is the same for P, D, and G, except for a rescaling of v with Γ. In particular, the curves reach the same minimal values and the minima are shifted to the right with increasing topology index Γ. Thus the structures P, D, and G perform equally well, but at different values of v. From Figure 6a, we see that all other structures cannot reach this low value of frustration since they have much broader distributions than P, D, and G (compare Table 3). Note that the scaling with v/Γ is valid in eq 12 for all structures and in eq 21 only for those related by a Bonnet transformation. Except for a factor of v/c02, the free-energy density (9) for r ) 0 is identical with 〈(Hl - c0)2〉ll2 (see eq 21). In Figure 7, we show the free-energy density of all considered phases for c0 ) 1/6, r ) 0, and κ/kBT ) 10. From the graph, we infer the phase sequence L-G-D-P with increasing water content; the other phases are not stable since they cannot reach that small an amount of bending frustration as can P, D, and G. As discussed above, eq 13 implies that the sequence G-D-P is determined by the topology index. Since the free-energy density for P rises again for small v, it would be more favorable for it to coexist with excess water. This mechanism is known as emulsification failure in surfactant systems31 and allows the sequence to be extended to include a P + W coexistence. The emulsification failure can be interpreted as a Maxwell construction of P with the pure water phase with f ) 0 at v ) 0. It follows from eq 9 and the scaling of 〈(Hl - c0)2〉ll2 with v/Γ that the free-energy density fb for P, D, and G can be written in the scaling form
fb ) vΩ(v/Γ)
(22)
with a universal function Ω(x). For small r, this function has Ω(0) ) 1, has a minimum at x = 4 c0/(1 + c0), and diverges at x = 4. A Maxwell construction shows that these structures always lie on a triple line, irrespective of the values for c0, r, and κ (compare inset of Figure 7). This means that if we calculated the whole phase diagram with the Maxwell construction, D and P would be stable only along a line. In the following we identify phase transitions between G, D, and P with the intersections of the free-energy curves, since under experimental conditions, small additional contributions to the free energy are certain to destroy the delicate free-energy balance resulting from the Bonnet transformation and then will lead to extended one-phase regions for D and P (or remove one or both of them from the overall phase diagram). In Figure 8 we use this reasoning to construct phase diagrams as a function of v and c0 for different values of r and κ/kBT. With increasing water concentration FW ) 1 - v, we always obtain the sequence L-G-D-P-W + P. The most stable IBCP is the gyroid structure G close to the lamellar phase at high values of v since it has both a narrow distribution of Gaussian curvature and the highest topology index. At lower concentrations of v, one finds narrow regions of stability for D and P since they have narrow distributions, too, but increasingly lower values of Γ. Increasing c0 shifts the loci of least frustration to higher v as evident from Figure 6a. Decreasing κj (i.e., increasing r) disfavors the IBCPs and increases the stability region for the lamellar phase, as can be seen by comparing parts a and b of Figure 8, where r is increased from 0.1 to 0.5, respectively. However, for c0 J r/2 there always will be a region of stability for the IBCPs, even for r ) 1. This stands in marked contrast to similar modeling for ternary systems where bicontinuous cubic phases are expected to disappear for r J 0.55.40 The explanation lies in the fact that changing the spontaneous curvature c0 of the monolayers amounts to changing the effective saddlesplay modulus of the bilayer; compare eq 4. Thus, increasing c0 can counteract the effect of making κj more negative. The bicontinuous phases are also favored by lowering κ without changing r, since this increases the fluctuation effects which favor IBCPs (positive renormalization of κj) and disfavor the lamellar phase (increase of
2094
Langmuir, Vol. 17, No. 7, 2001
Schwarz and Gompper
uniquely by the K distribution, which describes the internal geometry of a surface. This is evidenced by the self-intersections of the parallel surfaces for v J 0.8, which do not affect any of the considerations of the curvature energy described in the previous sections but would give a divergent contribution to the van der Waals energy. Thus, we conclude that the van der Waals interaction must break the Bonnet symmetry. We illustrate this point in Appendix A, where we approximate the TPMS by a face-center cubic (fcc) array of spherical shells. To estimate the frustration energy due to chain stretching, it has been suggested in ref 32 to calculate the variance of the distance between constant-mean-curvature companions of a given TPMS. For the D surface, the variance ∆L ) 〈(L - l)2〉 has been shown32 to be almost independent of v, with
∆L ) RDv2a2 ) 4RDA*2l2
(25)
and RD ) 0.000 35, where eq 6 has been employed in the second equality. Similar results can be expected for G and P, with prefactors RG and RP, respectively. We introduce a function R(A*), which is defined such that R(AP*) ) RP, R(AD*) ) RD, and R(AG*) ) RG for the three values of A* for the P, D, and G phases. Then, the stretching-energy density has the scaling form
Figure 8. Phase diagrams as a function of water volume fraction FW ) 1 - v and spontaneous curvature c0 for (a) r ) 0.1 and κ/kBT ) 10, (b) r ) 0.5 and κ/kBT ) 10, and (c) r ) 0.5 and κ/kBT ) 2. Two-phase coexistences are indicated by hatched regions. Part a is reproduced from ref 1. Copyright 2000 American Institute of Physics.
steric repulsion). This can be seen by comparing parts b and c of Figure 8, where κ/kBT is decreased from 10 to 2, respectively, while r is kept constant at 0.5. VI. Breaking the Bonnet Symmetry It is clear from our discussion in section II that all contributions to the free energy, which can be written in the scaling form vΩ(2l〈K〉1/2) ) vΩ(v/Γ), will preserve the triple-line coexistence of G, D, and P. Obviously, an interaction which introduces an additional length scale, in addition to l and 〈K〉-1/2, is needed to break the Bonnet symmetry and to lift the triple-line coexistence. Here we want to discuss van der Waals interactions and chain stretching. The nonretarded van der Waals interaction at distance d ) |r1 - r2| has the general form
w(d) ) -
AH 2
π
F1F2 d-6
(23)
where AH is the Hamacker constant and F1 and F2 are the number densities of atoms in the two interacting bodies. The total interaction energy between two colloidal particles of volumes V1 and V2 is then given by
W)
∫V
1
d3r1
∫V d3r2 w(|r1 - r2|) 2
(24)
Unfortunately, the energy density of the van der Waals interaction cannot be calculated easily in a cubic bicontinuous phase. However, it is clear from eq 24 that it is the embedding of the TPMS in three-dimensional space, which is important for the van der Waals energy. Distances between different parts of the bilayer are not determined
fs )
ks(l)l4 k (l)l4 l3 Es -2 2 s ) v ∆ l ) 2R(A*)A* v L 2κc02 V 2κc02 κc02 (26)
where Es is the stretching energy from eq 2 and eq 25 has been used for ∆L. Note that this energy density has been made dimensionless with the same factor as the curvature energy in section II; i.e., it is normalized to the free-energy density of the lamellar phase (compare eq 7). These results make the scaling argument of section II more precise. Equation 25 shows very clearly that the lattice constant a appears as a new length scale and that the Bonnet symmetry is broken in general by the stretching contributions. It is interesting to consider the possible functional dependencies of R(A*) in more detail. First, if R(A*)A*2 is a monotonically increasing (decreasing) function of A*, then P is stabilized (destabilized) by chain stretching due to its small value of A*, while D is disfavored (favored). In the case of increasing R(A*)A*2, D would not appear in the phase diagram, while in the case of decreasing R(A*)A*2, P cannot be present. Second, if R(A*)A*2 ) constant, the Bonnet symmetry is restored. Finally, R(A*)A*2 can of course be a nonmonotonic function of A*, which stabilizes the phase with the smallest value of R(A*)A*2. In any case, the difference between the three phases can be expected to be small, since the prefactor RD (and presumably also RG and RP) is very small; nevertheless, the stretching energy might account for small differences, such as the ones needed to produce the phase diagram of 2:1 LA/DLPC and water, which shows no triple-line symmetry but something close to it. A more precise statement requires the calculation of ∆L for G and P in the framework of the CMC-model. Other estimates of the stretching energy have been discussed in ref 21; they also indicate a breaking of the Bonnet symmetry. VII. Concluding Remarks In this work, we have investigated inverse bicontinuous cubic phases (IBCPs) in lipid-water mixtures. We systematically included all IBCPs which might be expected
Lipid-Water Mesophases
Figure 9. Experimental phase diagram for 2:1 lauric acid/ dilauroyl phosphatidylcholine and water as a function of water volume fraction FW and temperature T (adapted from ref 21). The temperature range from T ) 35 to T ) 50 °C corresponds roughly to the range 0.18-0.26 in spontaneous curvature c0 in the predicted phase diagram; compare Figure 8.
to be of physical relevance and found that only the ones with a narrow distribution of Gaussian curvature over the minimal midsurface are stable, i.e., P, D, and G. While this result has been anticipated in ref 39, we were able to prove it quantitatively for the first time by calculating the distributions of Gaussian curvature over the minimal midsurfaces of all relevant structures; from these we calculated the variance ∆ of the K distributions given in Table 3 which is the same for P, D, and Gsdue to Bonnet transformationssand higher for all other structures. The Bonnet transformation also implies that the freeenergy densities for P, D, and G have the scaling form fb ) vΩ(v/Γ) with hydrocarbon volume fraction v and topology index Γ. The topology index is a dimensionless inverse Euler characteristic which is independent of scaling and choice of unit cell; the higher its value, the less porous the structure. Since it is highest for the gyroid structure G, it has the largest region of stability. D and P have progressively lower values for Γ and have therefore smaller regions of stability at lower values for v. We also used the concept of the emulsification failure (which is well established in surfactant science) to explain the stability of IBCPs in excess water. A Maxwell construction with the excess phase shows that P, D, G, and excess water lie on a triple line for all values of c0 and r as a consequence of the Bonnet transformation. This “Bonnet symmetry” should not be interpreted as a degeneracy between G, D, and P, since their different topologies ensure that their regions of stability are arranged according to their respective values of the topology index. It should be noted that the free energy per unit area (rather than per unit volume) has the scaling form Ω(v/Γ); this might be interpreted as some kind of degeneracy (compare refs 21 and 39), since now all three structures attain the same minimal values. Our results agree nicely with experimental observations. At fixed temperature, it is always the sequence G-D-P which is observed with increasing water content, with G having the largest region of stability and only D or P coexisting with excess water.23 For example, for monoolein and water at 40 °C, with increasing water concentration the sequence L-G-D has been found, where G has the largest region of stability and D coexists with excess water.23 Our phase diagram for r ) 0.1 and κ/kBT ) 10 (Figure 8a) agrees particularly well with the one for 2:1 LA/DLPC and water21 (compare Figure 9). We do not treat the inverse hexagonal phase, which experimentally occurs at large v and large c0, but we showed above that IBCPs will occur at large c0 and small v, and argued that at large v stretching energy should become comparable to bending energy. This, in fact, favors the inverse hexagonal phase, whose geometry requires chain stretching.
Langmuir, Vol. 17, No. 7, 2001 2095
Finally, we can discuss the role of temperature for 2:1 LA/DLPC and water by using the ansatz c0 ) b(T - Tb), with a balanced temperature Tb and a phenomenological constant b. Using l ) 11.4 Å, we find b ≈ 0.0011/ÅK and Tb ≈ 0 °C. Since the chain melting temperature is known to be Tm ) 30 °C, we conclude that the balanced temperature Tb is well below the main transition. This seems to be the generic case for lipid systems; it might be related, e.g., to the anomalous swelling of PC bilayers when the main transition is approached from above.57 It is interesting to note that the constant b obtained here for 2:1 LA/DLPC and water is very close to the corresponding value b ) 0.0012/ÅK for monolayer of C12E5 in ternary microemulsions.46 In this paper, we only considered inverse bicontinuous phases with cubic symmetry. A similar analysis for noncubic phases should be straightforward, since here Weierstrass representations are known for the main structures, especially for the rhombohedral and tetragonal variants of G, D, and P.27 However, we expect that these phases will not be more favorable than G, D, and P for the same reason as for the other cubic phases. To make further progress in understanding their phase behavior, it seems more important to investigate the different mechanisms which might break the Bonnet symmetry between G, D, and P. In our view, two calculations are needed now: (i) of the van der Waals energy of G, D, and P and (ii) of the variance of the distance distribution of a constant-meancurvature model of G and P. Appendix A: Estimate of van der Waals Energy Density For a rough estimate of the energy density due to the van der Waals interaction in cubic bicontinuous phase, we use a fcc array of spherical shells of radius R, thickness 2l, and lattice constant a. Here, the lattice constant a corresponds to the conventional unit cell and is taken to be the same as for the corresponding bicontinuous phase, and R is determined by the requirement of equal area A per unit cell, so that A* ) 16 π(R/a)2. Thus we have a/l ) 2A*/v again. The distance of closest approach between two neighboring spheres is D ) 21/2a/2 - 2R; in order to avoid contact between them, D > 0 is necessary, which implies A* < 2π; this inequality is satisfied for all TPMS studied here (compare Table 1). For the crystal, we consider the interaction between nearest-neighbor spheres only. The van der Waals energy between two hollow spheres is described by the Girifalco potential, which is often used to model the van der Waals interaction between two buckyballs C60.58 For our purpose it is sufficient to use the Derjaguin approximation for the Girifalco potential, which describes the leading contributions for close approach correctly and shows that the van der Waals attraction for two hollow spheres scales as W ∼ AHRl2/D3. We thus find for the van der Waals energy density
fvdW )
() (
AH l l3 W ∼ 2 V 2κc0 κc02 a
5
)
R 1 R a 21/2 a
-3
(A1)
This energy density has again been made dimensionless with the same factor as the curvature energy in section II. Its scaling shows that two ratios of length scales, (56) Templer, R. H. Curr. Opin. Colloid Interface Sci. 1998, 3, 255263. (57) Richter, F.; Finegold, L.; Rapp, G. Phys. Rev. E 1999, 59, 34833491. (58) Girifalco, L. A. J. Phys. Chem. 1992, 96, 858-861.
2096
Langmuir, Vol. 17, No. 7, 2001
Schwarz and Gompper
R/a ∼ (A*)1/2 and a/l ∼ A*/v, determine the van der Waals energy. It is intersting to note that the curvature-energy density in this model of spheres scales as
[
fb ) v 1 -
1 c0(R/l)
]
2
Since R/l ∼ (A*)3/2/v ∼ Γ/v (for fixed χ), eq A2 is consistent with the scaling function (22). This indicates that our identification of the sphere radius and lattice constant gives a consistent approximation for both van der Waals and bending energies.
(A2) LA0013805
