Originally published as Biophys J. BioFAST on March 13, 2006.
doi:10.1529/biophysj.105.062778
Biophysical Journal 90:4104-4118 (2006)
© 2006 The Biophysical Society
Effect of Hydrophobic Mismatch on Phase Behavior of Lipid Membranes
Elizabeth J. Wallace * 
,
Nigel M. Hooper * and
Peter D. Olmsted
* School of Biochemistry & Microbiology, Institute of Molecular and Cellular Biology and Leeds Institute of Genetics, Health and Therapeutics, School of Physics & Astronomy, University of Leeds, Leeds, United Kingdom
Correspondence: Address reprint requests to P. D. Olmsted, School of Physics & Astronomy, University of Leeds, Leeds, LS2 9JT, United Kingdom. Tel.: 44-113-343-3830; E-mail: p.d.olmsted{at}leeds.ac.uk.
 |
ABSTRACT
|
|---|
We investigate the competing effects of hydrophobic mismatch and chain stretching on the morphology and evolution of domains in lipid membranes via Monte Carlo techniques. We model the membrane as a binary mixture of particles that differ in their preferred lengths, with the shorter particles mimicking unsaturated nonraft lipids and the longer particles mimicking saturated raft lipids. We find that phase separation can be induced upon increasing either the ratio 
of the hydrophobic surface tension J to the compressibility modulus 
. J/ determines the decay length for thickness changes. When this decay length is larger than the system size the membrane remains mixed. Furthermore, increasing the thickness relaxation time can induce transient phase separation.
|
INTRODUCTION
|
|---|
Strong evidence has been obtained that suggests the presence of heterogeneities in the plasma membrane of many cells (1
–5). One type of membrane heterogeneity, termed rafts, is enriched in cholesterol, sphingomyelin (SM), and certain membrane proteins (6,7). Rafts have putative roles in many physiological processes, such as signal transduction, endocytosis, apoptosis, protein trafficking, and lipid regulation (7–11).
Raft lipids typically have saturated hydrocarbon chains. For example, SM comprises a sphingoid base that has a saturated hydrocarbon chain plus a very long amide-linked saturated acyl chain that is, on average, 20–24 carbons in length (12). Glycerol-based phosphatidylcholine, on the other hand, which is the major class of nonraft lipids that have, on average, 18 carbons per acyl chain (13), is less saturated than naturally occurring sphingomyelins, which generally contain one saturated and one unsaturated acyl chain (14). These unsaturated bonds produce kinks in the lipid chains, increasing the area per molecule. Therefore, this leads to a reduction in lipid thickness due to volumetric considerations. Cholesterol, a molecule that is enriched in lipid rafts, has a shorter hydrophobic length of 17.5 Å (15). However, it does not change the length of SM. With or without cholesterol, C18:0 SM has been found to form bilayers with a thickness of 46–47 Å (16). Hence, because of the differences in lipid thickness between raft and nonraft lipids, and since energetically it can be expected that the lengths of the hydrophobic moieties of neighboring membrane components will be approximately equal to avoid unfavorable exposure of hydrophobic surfaces to a hydrophilic environment, it is reasonable to assume that nonraft lipids should coalesce and predominantly give rise to smaller bilayer thicknesses compared to raft lipids (17,18). Indeed, in studies of model membranes comprising dioleoylphosphatidylcholine (DOPC), SM, and cholesterol, phase separation is observed giving rise to SM- and cholesterol-rich raft domains (19–21). Atomic force microscopy studies reveals that these domains are typically of order 1-nm thicker that the surrounding DOPC-rich region (19).
However, it has been shown that proteins are able to integrate into a membrane while tolerating a length variation of 
10 amino acids, depending on the amino acid composition (22,23). This therefore raises the questions: How important is hydrophobic mismatch for membrane structure and organization, and how do membranes relieve hydrophobic mismatch between proteins and lipid bilayers? Possible adaptations for the relief of hydrophobic mismatch in membranes include phase separation of lipids of different thicknesses, ordering and disordering of the acyl chains, peptide backbone deformation, peptide aggregation, peptide tilt, no transmembrane association of the peptide with the membrane, or nonlamellar phase formation. Little is known about which will be favored under certain conditions (24,25).
The goal of this article is to investigate how hydrophobic mismatch in a simple model membrane comprising two species of particle influences the membrane phase behavior. Specifically, we model adaptations to hydrophobic mismatch by stretching or compressing the particles, hence ordering or disordering the lipid acyl chains.
Theories that address acyl-chain stretching or compression due to the incorporation of a peptide of a different hydrophobic thickness include the phenomenological approaches by Owicki et al. (26,27) and Jähnig et al. (28–30). In these theories, they assume the location of the proteins to be fixed and so they do not induce thermodynamic phase separation. The phenomenological mattress model by Mouritsen and Bloom (31) is a two-component real solution theory and hence allows for phase separation. In their model, they relate the energy stored in the undulations of the membrane surface caused by the mismatch to the elastic properties of the lipids and proteins. They do not include microscopic detail of the lipids, but use as input the known thermodynamic properties of the pure lipid system. They also include indirect lipid-protein interactions induced by the mismatch as well as direct lipid-protein van der Waals-like interactions between the hydrophobic parts of the lipid bilayer and the proteins. The mattress model has been replicated in a Monte Carlo simulation scheme by Sperotto and Mouritsen (32). They allowed for different microstates of the lipids, classified according to Pink's 10-state model (33), hence enabling a pure lipid bilayer phase transition. A major theoretical advance was the work of Fattal and Ben-Shaul (34), who provided a molecular theory for the behavior of the lipid chains. This molecular modeling was combined with phenomenological free energy contributions accounting for the opposing effects of headgroup repulsions and hydrocarbon-water surface tension. Duque et al. (35) described the effects of an embedded protein in a bilayer via molecular theory, which yielded the free energy of the entire system. All these theories predict a significant stretch of the acyl chains of lipids that are adjacent to a long peptide. These perturbations in lipid thickness decrease for the lipids that are at a greater distance from the peptide.
Many experimental studies of model membranes support these theories. Deuterium (2H) NMR measurements on lipids with perdeuterated acyl chains show that a reduction in both acyl-chain order and bilayer thickness occurs when an embedded peptide has a smaller hydrophobic thickness than that of the lipid bilayer (36,37). Additionally, electron cryomicroscopy measurements showed a reduction in bilayer thickness when a short peptide was incorporated into a DOPC liposome (38). For the opposite case when a peptide has a hydrophobic thickness that is greater than that of the bilayer, it has been found via (2H) NMR experiments that the bilayer increases its acyl-chain order and thickness. Interestingly, upon increasing the peptide hydrophobic length, the membrane thickness does not match this increment, although its thickness does always increase (36,37,39). Hence, it appears that complementary mechanisms of adaptation to hydrophobic mismatch occur concomitantly with changes in lipid acyl-chain order and thickness.
Issues that have not been addressed previously relating to acyl-chain stretching or compressing to relieve hydrophobic mismatch include the initial instability of a membrane that could, for example, be induced biologically after delivery of lipids having a different hydrophobic thickness to a membrane, and also the potential induced large-scale structure. In this article, we study how the hydrophobic mismatch between two species of particle in a bilayer influences the membrane phase behavior and the kinetics of phase separation. In particular, a large hydrophobic mismatch between particles that are stiff can induce phase separation since the particles are unable to stretch or compress to relieve the mismatch. The two species of particle, U and S, represent unsaturated and saturated lipids, respectively, with the U particles having a shorter preferred hydrophobic length than the S particles. The mechanism of adaptation to hydrophobic mismatch is modeled by stretching and compression of the particles, analogous to stretching and compression of lipid acyl chains. To implicitly treat a bilayer, a particle is assumed to be the same lipid on opposite sides of the bilayer. Hence, a change in particle hydrophobic thickness will stretch (or compress) both lipids in opposite directions. We investigate the effects of hydrophobic mismatch on membrane phase behavior in two ways. Firstly, we present a mesoscopic calculation that probes the intermediate to late stages of composition and bilayer thickness growth. These regimes are explored via Monte Carlo computer simulation. Secondly, we present an analytic treatment describing the early stages of growth.
|
MESOSCOPIC MODEL
|
|---|
Free energies
Initially we study the effects of hydrophobic mismatch on domain morphology for a mesoscopic model. We model the bilayer as a two-dimensional square lattice with sides of length L. Each lattice site, having an area a2, is occupied by either a U or S particle. The particles have two degrees of freedom. They can laterally exchange positions with a neighboring particle or they can change their thickness. A particle is assumed to be the same lipid on opposite sides of the bilayer and so a change in particle thickness will stretch (or compress) both lipids in opposite directions. All energies are measured in units of kBT and all lengths are measured in units of the lattice spacing a.
The free energy G of the membrane comprises three terms: the lipid-lipid interaction energy Gint, the lipid stretching energy Gstretch, and the hydrophobic mismatch energy Gmismatch.
The lipid-lipid interaction energy is given by
 | (1) |
where 
i
is equal to 1 if species is at lattice site i or equal to 0 if species is not at lattice site i. The value can be U or S. The value Vß is the contact energy of neighboring species and ß. The physical contribution to Vß is from electrostatic and van der Waals interactions between the lipids. However, these interactions are not modeled explicitly. This term can lead to phase separation in a two-component system if the strength of the energetic interaction between U and S particles (VUS) relative to their self-interactions (VUU, VSS), 
, given by
 | (2) |
satisfies > MF = 2 within mean-field theory (40), or > c = 3.526 in a physical system incorporating critical fluctuations (Ising model) (41). Hence, for VUU = VSS = 1, phase separation occurs if VUS > 1.88.
The lipid stretching energy is given by
 | (3) |
where li is the actual thickness at site i, and l0 is the preferred thickness of species at site i. The preferred length of a particle is constant throughout a simulation. Hence, phase transitions between liquid-ordered and gel states are not modeled since they are usually accompanied by changes in lipid length. The value 
is the compressibility modulus of particle . The stretching energy arises due to a lipid not having its preferred thickness. For simplicity 
and 
are assumed to be equal in this study, hence we will refer to as . However, in a physical system, will be larger for lipids that are more saturated, therefore having straighter and less flexible hydrocarbon chains.
The hydrophobic mismatch energy is proportional to the exposed hydrophobic area and is given by
 | (4) |
where the exposed hydrophobic area is linear in the particle thickness difference. We have not included an explicit bending penalty in the mesoscopic model because bending, a long wavelength phenomenon, will arise from coarse-graining Eqs. 3 and 4 up to longer length scales. The value Jß is the hydrophobic surface tension between nearest and next-nearest species and ß. The value Gmismatch will be minimized when the particles are of equal thickness. Physically, this means that the lipid headgroups will be directly adjacent, hence avoiding contact between the aqueous environment and the lipid hydrocarbon chains, and/or avoiding contact between the water molecules around the lipid headgroup of one lipid and the hydrocarbon chains of the other lipid. In this study, we do not model lipid headgroups. Hence, we assume that the hydrophobic surface tension between a lipid headgroup and hydrocarbon chains, and the hydrophobic surface tension between the aqueous environment and the hydrocarbon chains, are equal. Additionally, for simplicity we assume Jß 
J independent of particle type. However, in a physical system the hydrophobic surface tension will be higher for unsaturated lipids due to the acyl chains being kinked. This leads to a larger available contact area between the lipid chains and water around neighboring lipid headgroups or the aqueous environment.
Model approximations
Our model assumes that the two monolayers are symmetric. Hence, a concentration of a short component in one leaflet implies a similar concentration in the opposite leaflet (Case I), and thus local membrane thinning. Alternatively, an excess of short components in one leaflet could be compensated by an excess of longer components in the opposite leaflet (Case II), thereby reducing the hydrophobic mismatch. There are two main contributions to the difference in free energies 
G = GI – GII between phase-separated morphologies in these two cases. The first contribution is due to the hydrophobic mismatch between domains. This incurs a free energy penalty Gline 
Ld, where Ld is the domain interface length and is the line tension. The second contribution is due to differences in the free energy of bilayer assembly 
, since we expect the symmetric bilayer to have the lower energy state. The competition between these two free energy differences leads to a length scale, 
. At early times, 
, suggesting that Case II will be the preferred configuration. However, at later times Ld will exceed 
, which should favor Case I. Our model is appropriate for the limit of 
. We leave the very interesting physics of the interplay between these two configurations to future work.
We also assume negligible frame tension, which prevents the area per particle from decreasing. The work done against a frame tension 
is given by
 | (5) |
where v is the volume of species . This would therefore shift the phase boundaries slightly; i.e., a deeper quench will be needed to effect phase separation that would lead to an overall smaller area per particle.
Simulation details
The simulations are initialized by distributing the particles randomly, with each species assigned its preferred thickness. We only consider the symmetric case of equal compositions U = S = 1/2 here. The thicknesses then undergo a preliminary relaxation period, where each particle is allowed to relax a possible 1000 times. A thickness move consists of randomly selecting a particle and choosing a thickness change randomly between ±
lmax; this change is accepted or rejected according to the Metropolis criterion (42). The choice lmax = 1 allowed thermal fluctuations to be sampled efficiently. After the initial thickness relaxation, composition moves are made in conjunction with thickness relaxation. Lateral exchanges are chosen that preserve the thicknesses of the individual particles: e.g., a U particle of thickness l at site i and an S particle of thickness l' at site j are swapped so that site i now contains an S particle of thickness l'. These exchanges are implemented by the usual Kawasaki dynamics (43) and accepted according to the Metropolis criterion. After each possible composition exchange, the thicknesses on sites i, j, and those of the particles in the surrounding shell of lattice sites, are allowed to relax a possible nr times per particle. The thickness relaxation is performed regardless of whether the particle lateral exchange was accepted or rejected, to allow continual thickness relaxation. We define one Monte Carlo (MC) cycle as the number of steps required for each lattice site to have the opportunity to have a particle lateral exchange. The ratio of characteristic times for thickness relaxation and diffusion can be adjusted by changing nr, with a decrease in nr leading to slower thickness relaxation.
|
MORPHOLOGY AND EVOLUTION OF DOMAINS USING THE MESOSCOPIC MODEL
|
|---|
Domain morphologies
Simulations were performed to investigate the effects of hydrophobicity on domain morphology. Firstly, however, we will demonstrate the effects of in the absence of coupling to thickness for comparison purposes. Fig. 1 a shows the domain morphologies obtained after 131,072 MC cycles in the absence of coupling composition with thickness (
) and increasing from 0 to 8. For = 0 the membrane remains mixed since there is no driving force for phase separation. As is increased, the system moves toward the critical point between the one-phase and two-phase regimes. Hence, larger scale fluctuations occur. As the system crosses the critical point (c = 3.526), phase separation occurs. Upon further increasing , purer phase-separated domains are observed that coarsen more slowly since the energy cost of a U–S contact becomes higher.
When considering hydrophobicity we find two ratios of interest, and 
, with having units of length while has units of energy. To compare the effects of hydrophobicity with the effects due to (Fig. 1 a), we use the ratio since this has units of energy, similarly to Vß (recall that 
). We anticipate similar trends to those observed in Fig. 1 a, provided the particles U and S have different preferred thicknesses. For small , the membrane should remain mixed since the energetic cost of a difference in thickness between neighboring particles will be low and hence there will be no driving force for phase separation. Increasing should drive the system into the phase coexistence regime to minimize the energy of the system by placing particles with similar thicknesses together.
As can be observed in Fig. 1 b, the effects of hydrophobicity probed by increasing are indeed similar to the effects of increasing . Domains are not observed for 
, since there is no impetus for phase separation. As is increased to 0.0875, critical fluctuations are observed. For 
, phase separation occurs, where domains coarsen more slowly upon increasing , due to higher incompatibilities between the unlike species.
We next investigate the effects of increasing the length for = 0 and 
for a system of size L = 50 (Fig. 1, c and d). As can be observed, the phase behavior differs significantly at high to that observed at high . Domain formation does not occur for 
and = 0, since there is no incentive for phase separation. As is increased to 0.2, critical fluctuations are observed. Upon increasing from 0.4 to 2, phase separation occurs. A similar trend is observed to increasing either (Fig. 1 a) or (Fig. 1 b) in that as the system moves further from the critical point, domains coarsen more slowly due to higher incompatibilities between the unlike species. However, as is increased from 2 to 53, domains become larger and less pure after the same number of MC cycles. Hence, there is a reversal in behavior. In the next subsection, we will discuss the relative rates of growth in composition and thickness, and their dependence on the relative relaxation times for composition and thickness. Interestingly, for 
, the membrane mixes for all initial configurations. To understand this behavior, simulations were performed for phase-separated membranes where only the particle thicknesses were allowed to relax (Fig. 2), therefore probing whether phase separation will be obtained for a given value of at late times. For and L = 50 (Fig. 2 c), the thickness profile is approximately flat after equilibration, implying that the composition will decouple from the thickness for these simulation parameters, which for = 0 would imply a mixed state. Indeed, as is increased to 100 in Fig. 1, c and d, the particle thicknesses become equal and, hence, there is no impetus for phase separation from hydrophobicity. However, upon increasing the system size L to 100 for (Fig. 2 d), there is a clear gradient in the thickness profile at the domain interface after equilibration. The difference in thickness profiles for different system sizes suggests that there is a decay length 
, due to the competition between hydrophobic mismatch and particle stretching, which controls phase separation. For 
decay > L, hydrophobic mismatch dominates, resulting in particle thicknesses becoming equal and therefore decoupling from the composition. However, for decay < L, particle stretching dominates, which can result in phase separation depending on the value of . For and L = 50, decay > L and therefore phase separation is not observed. Note that starting the simulation from a phase-separated and flat state for this set of parameters results in remixing. Hence, the lack of phase separation at high is not thought to be a kinetic effect, but a finite size effect.
View larger version (14K):
[in this window]
[in a new window]
|
FIGURE 2 Thickness profiles taken at L/2 for simulations that had a fixed phase-separated composition allowing the particle thicknesses to fully equilibrate. Each phase-separated domain was composed of either all U or all S, and had sides of length L and L/2. , lu0 = 36, ls0 = 54, lmax = 1.0, u = s, and = 0. (a)  and L = 50, (b) and L = 100, (c) and L = 50, and (d) and L = 100. Time t is in MC cycles, where in this figure an MC cycle is defined as the number of steps required for each lattice site to have the opportunity to have a particle-thickness change.
|
|
We next determine the versus phase diagram (Fig. 3) for = 0 and system size L = 50. As can be seen, as approaches the system size L, finite size effects lead to a significant increase in the magnitude of required to induce phase separation.
Fig. 1 d shows the membrane morphologies obtained when increasing for = 5 and . For the membrane has phase-separated due to . Upon increasing to 0.4, domains decrease rather than increase in size. This is due to the membrane already being below the critical point at . Above 
, similar trends are observed to those shown in Fig. 1 c, with the composition decoupling from thickness at since decay > L.
Fig. 4 shows the versus phase diagram for symmetric quenches (i.e., u = s) into the phase coexistence regime for L = 50 and for (note that, as with Fig. 1, c and d, the value of is not constant). Above c = 3.526, the membrane will demix regardless of the value of . However, below c, the membrane will mix for either small or large values of . The lower limit of that leads to phase separation decreases upon increasing , and the upper limit of that leads to phase separation increases upon increasing . As discussed above, this upper phase boundary between mixing and phase separation for < c is a finite size effect due to decay > L. For = 0, increasing at fixed in Fig. 4 traces a parabola in Fig. 3, which intersects the phase boundary in Fig. 3 at two points. This leads to the progression of states from one-phase to two-phase to one-phase.
View larger version (6K):
[in this window]
[in a new window]
|
FIGURE 4 versus phase diagram for , lu0 = 36, ls0 = 54, lmax = 1.0, nr = 1, L = 50, u = s, and simulation time t = 131,072 MC cycles. The crosses (x) are simulation data marking the approximate boundaries of the two-phase regime. The dashed-dotted line marks the upper phase boundary between the one- and two-phase regimes, where the position of this line depends on the system size. The dotted line is predicted behavior and is only schematic. For > c = 3.526 (dashed line), phase separation will occur regardless of the value of J. For comparison purposes between the phase diagram shown here and the phase diagram in Fig. 3, the lower phase boundary for = 0 occurs at  , while the upper phase boundary occurs at  .
|
|
The results presented thus far illustrate trends in behavior upon changing and the hydrophobicity parameters J and . We will now investigate membrane phase behavior for biologically relevant parameters. X-ray diffraction measurements of lipid bilayers have determined that at 30°C, fluid phase DOPC has a cross-sectional area of 73 Å2 and a bilayer thickness of 37 Å (44). At 22°C, gel phase SM (C18:0) has been found to have a cross-sectional area of 45 Å2 and a bilayer thickness of 50.5 Å (45). Assuming no area change upon mixing, 1:1 SM/DOPC bilayer has an average cross-sectional area of 59 Å2, giving rise to an average distance between neighboring lipids of 8.6 Å. Lengths in the simulations are measured in units of the lattice spacing a. Therefore, reasonable estimates for the preferred thickness of the U species of particle lu0 that mimics DOPC and the preferred thickness of the S species of particle ls0 that mimics SM are 4.3 a and 5.9 a, respectively.
Measured values for the area compressibility modulus 
A on fluid lipid bilayers and free biological cell membranes range from (100–230) mJm–2 = (0.242–0.556) kBT Å–2 at 300 K (46–48). The elastic energy per molecule can be expressed as a difference between the actual area A and the preferred area A0 per molecule by (49)
 | (6) |
and in the mesoscopic model studied here, the elastic stretching energy per molecule is expressed as a difference in thicknesses by
 | (7) |
If the change in volume v = Al of the lipid as it stretches is negligible, then
 | (8) |
Therefore, an average value of for DOPC and SM is 0.01 kBT Å–2 = 0.7 kBT a–2.
The hydrophobic surface tension J, due to tail-water or tail-headgroup interactions, provides the dominant contribution to the energy cost due to a mismatch in thicknesses of neighboring lipids. The interfacial tension 
WC of hydrocarbons with bulk water lies in the range (40–50) mJm–2 = (0.097 – 0.121) kBT Å–2 at 300 K (49). This tension will mainly be due to the hydrophobic effect. When the headgroup is in contact with the lipid tails, WC needs to be scaled by the amount of water present in the heads, which, for a DPPC bilayer, is approximately one-third of the total headgroup volume (44). Hence, a rough estimate for the interfacial tension between a lipid headgroup and lipid acyl chains HC is 0.04 kBT Å–2. However, entropy that is associated with the hydrogen-bond network of water contributes on the order of 85% to the hydrophobic interaction (49). Since the hydrogen-bond network of water will be significantly perturbed in the lipid head environment, the interfacial tension should be even lower, of the order of 85%. This leads to HC 0.006 kBT Å–2. In the simulations, J is given by
 | (9) |
Therefore, for an average cross-sectional area of DOPC and SM, J = 0.05 kBT Å–1 = 0.43 kBT a–1, assuming that hydrophobic mismatch only occurs between lipid heads and neighboring lipid tails.
Fig. 5 illustrates the effects observed when simulating a membrane comprising two species of particle mimicking SM and DOPC, with the rough estimates for J and . As can be seen, phase separation does not occur during the timescale of the simulation; however, increasing by a factor of 10 results in domain formation.
View larger version (42K):
[in this window]
[in a new window]
|
FIGURE 5 Domain morphologies (a) and thickness profiles (b) taken at L/2 for simulations having approximate biologically relevant parameters with lu0 = 4.3, ls0 = 5.9, lmax = 1.0, nr = 1, L = 50, u = s, t = 131,072 MC cycles, and = 0 (rough estimates for biologically relevant and J calculated in Morphology and Evolution of Domains Using the Mesoscopic Model are 0.7 and 0.43, respectively). U (solid); S (shaded). The x and y axes range from position x = 0...50 and thickness l = 4...7, respectively, for each graph.
|
|
Domain growth and effect of changing relative relaxation times between thickness and composition
We will next examine the growth in composition and thickness. Fig. 6 shows both the length scale of compositional growth * (determined from the peak in the structure factor) and the growth in the root mean-square thickness fluctuations lRMS for various values of J for nr = 1, u = s, and . Increasing J from 0.04 to 0.2 leads to slower coarsening with a change in power law, indicative of a change in the mechanism of growth. Interestingly, initial domain growth is not observed for J equal to 0.8, 3.5, and 5.3, with the delay in growth increasing for larger values of J. This is due to increasingly slower relaxation of the particle thicknesses. Note that since all quenches performed were symmetric, spinodal decomposition will have been the mechanism of domain growth. Growth during the early stages of such domain coarsening is expected to be due to diffusion, while coarsening at the late stages is expected to be due to interface-driven hydrodynamics (50). In this study, only the diffusion is simulated. Therefore, our model data will most accurately reproduce domain growth in atomic force microscopy experiments on supported lipid bilayers, in which growth by hydrodynamics is expected to be significantly damped. This is because the thin hydration layer between the support and the lower leaflet of the bilayer will exert strong hydrodynamic stress.
View larger version (14K):
[in this window]
[in a new window]
|
FIGURE 6 Characteristic length scale of growth in composition * (a,b) and growth in RMS thickness fluctuations lRMS (c,d) for J from 0 to 10 (for J from 0.0 to 0.2, each trace is an average of 10 data sets; and for J from 0.8 to 10.0, each trace is an average of five data sets) for , lu0 = 36, ls0 = 54, lmax = 1.0, nr = 1, L = 50, u = s, and = 0. There is no compositional growth for J = 0 and 10.
|
|
The dependence of the delay of compositional growth on thickness relaxation for J = 0.8, 3.5, and 5.3 can be demonstrated by changing the relative relaxation times between composition and thickness. This can be done by changing nr. For example, an increase in nr leads to a faster thickness response, or alternatively a relative decrease in particle lateral mobility. Fig. 7 illustrates the effects of changing nr on both the growth in composition and thickness, where each growth trace is averaged from five simulations. For small values of J (J = 0.04), increasing nr from 1 to 4 has little effect on growth in composition and thickness since the impetus for equilibrating the particle thicknesses is also small. Increasing J to 3.5 leads to faster thickness relaxation for nr = 4 than for nr = 1. This increase in thickness relaxation enables the growth in composition to also hasten. Again, for J = 5.3 the thickness response is faster for nr = 4 than for nr = 1. However, in this case, there is no overall compositional growth observed for nr = 4 during the timescale of the simulations.
View larger version (10K):
[in this window]
[in a new window]
|
FIGURE 7 Effects of changing the relative thickness relaxation time on the characteristic length scale of composition growth * (a–c) and RMS thickness fluctuation growth lRMS (d–g) for J equal to 0.04, 3.5, 5.3, and 10.0 (each trace is an average of five data sets). , lu0 = 36, ls0 = 54, lmax = 1.0, L = 50, u = s, and = 0. There was no compositional growth for J = 5.3 and nr = 4, and J = 10.0, nr = 1, and nr = 4.
|
|
The magnitude of determines whether or not an increase in nr enhances compositional growth. For smaller , the stretching free energy Gstretch will dominate the free energy penalty, which will encourage particles to have their preferred thicknesses. This occurs faster for larger nr. Therefore, within a concentration fluctuation that enhances a given species, the particle thicknesses will quickly relax toward their preferred thicknesses for large nr, minimizing both the stretching free energy (in the bulk) and also the hydrophobic mismatch free energy (which only occurs at the interface). Hence, an increase in nr will lead to faster compositional growth for small . However, for large the hydrophobic mismatch free energy Gmismatch dominates the free energy penalty. Gmismatch will encourage the particle thicknesses to become equal regardless of composition, with this occurring faster for larger nr. Hence, because the thickness gradient will be smaller for larger nr, thermal fluctuations will be able to break up concentration fluctuations that lead to clustering of like species more easily. Indeed, in Fig. 8 (J = 5.3 and nr = 4), the composition profile taken at t = 65,536 MC cycles clearly shows the early stages of phase separation, while at a later time (t = 131,072 MC cycles) the membrane has remixed. It should be noted that both the growth in composition * and the growth in thickness fluctuations lRMS shown in Fig. 7 are averaged from five data sets. Therefore, the compositional growth for J = 5.3 (Fig. 7) does not show this transient domain growth. It is anticipated that for smaller nr, these transient domains will have a longer lifetime since the thickness gradient at the domain interface will take longer to relax. In a physical system it is expected that the thickness relaxation will be fast, implying short-lived transient domains.
Thus far we have demonstrated via Monte Carlo simulation that hydrophobic mismatch between neighboring particles can lead to phase separation. Two ratios of interest have been identified; the ratio of the hydrophobic surface tension J to the compressibility modulus , which has units of length, and the ratio , which has units of energy. Increasing either of these ratios leads to phase separation, with the phase-separated domains coarsening more slowly upon increasing the ratios above zero. However, upon further increasing , domain sizes become larger after the same number of MC cycles. Interestingly, for high enough the membrane remains mixed; this is a finite size effect, due to a decay length that becomes greater than the system size L at high . Physically this is because the gradient in the particle thicknesses at the domain interface is too small to maintain phase separation after a concentration fluctuation. It has also been shown that faster thickness relaxation affects the speed of compositional growth.
|
CONTINUUM MODEL
|
|---|
In the previous section we examined coarsening and the late stages of phase separation due to hydrophobic mismatch. We will now use a continuum model to probe the initial instability of a membrane due to hydrophobic mismatch and the resultant initial domain growth.
The continuum model describing the evolution of composition and lipid thickness l can be written, to fourth order in the gradients of composition and thickness, as
 | (10) |
where the term containing 
is the coarse-grained version of Eq. 3. The value is related to the compressibility modulus in the mesoscopic model by 
. The value l0() is the preferred local lipid thickness. We expect coarse-graining to lead to a l0() that depends on 
. For large , l0() is expected to be a very weak function of . The value f0() is the free energy per area of mixing, given by
 | (11) |
The terms involving gradients in Eq. 10 are coarse-grained versions of Eqs. 3 and 4. Specifically, the terms comprising gl, gl, and gll are coarse-grained versions of the hydrophobic mismatch (Eq. 4), with the terms containing gll and gl describing the bending penalty (51). The values g and gl are the second-order gradient coefficients in composition and thickness, respectively, and gl, g, and gll are the next terms in the gradient expansions for composition and thickness, respectively. The higher order g and gll gradient terms are required to stabilize the growth rate at small length scales.
Phase separation
The local composition variable (typically area fraction) obeys the continuity equation (mass conservation),
 | (12) |
where the flux 
of material is given by Fick's law (40),
 | (13) |
M is the particle lateral mobility and µ G/ is the chemical potential. Upon combining Eqs. 12 and 13 we can expand the composition dynamics to linear order in deviations of composition and thickness. We use a Fourier expansion,
 | (14) |
 | (15) |
and we assume an initial configuration of uniform composition 0 and a corresponding average thickness l0.
To linear order in the q 
0 modes, we find
 | (16) |
where 
is a small perturbation in composition, 
, and 
.
Thickness growth
The thickness evolves by relaxational kinetics,
 | (17) |
where the friction coefficient 
is due to dissipation within the membrane. Using Eqs. 10 and 17, we find
 | (18) |
In Fourier space this is given by
 | (19) |
Coupled composition and thickness growth
Finally, we consider simultaneous composition and thickness evolution after a quench. Equations 16 and 19 can be written as
 | (20) |
where
 | (21) |
The value 
controls the dynamic coupling. For < 1, the diffusive coarsening is faster than thickness growth, while for > 1 the thickness evolves faster than diffusive coarsening. Note that > 1 is likely to be the more physical regime since the frictional forces incurred upon lipid stretching or compressing will be less than those incurred upon lipids exchanging lateral positions. > 1 is indeed found in the following simple calculation: For lipids in a bilayer the diffusion coefficient, D 5 x 10–8 cm2 s–1 (52), is roughly estimated as D 
a2/2
d, where a is a molecular diameter and d is the hopping time. Hence, for a 9 Å we estimate d 80 ns. The relaxation time l for the slowest peristaltic mode of a lipid bilayer, which corresponds to the relaxation of thickness fluctuations, calculated via molecular dynamics simulation is 4 ns (53). Hence, the ratio of these times d/l is 20.
Assuming the solutions
 | (22) |
 | (23) |
we have
 | (24) |
The eigenvalues (
1, 2) of the matrix V yield the rates 
q, which govern the growth of composition and thickness fluctuations. Because the second eigenvalue 2 was stable in the entire q-range for the parameters we consider, we will focus on the effect of 1 on domain growth.
Parameters are chosen as follows. All lengths, including the thickness variable l, are scaled by the lattice size a, and all energies are scaled by 
. The gradient terms g, gl, and gll can be, in principle, derived as higher-order expansions of terms 
, etc. Hence we write these as
 | (25) |
where and l are the ranges of the interactions in the gradient expansion. Since all gradient expansions will be governed by the lattice size a, which is of the order of a lipid diameter, we choose = l = a 9 Å. Note that the gradient expansion in thickness may also involve the thickness l0, but the largest effect should be due to splay near the surface, so we take as an estimate l = a (an upper bound on this may be 
). This leads to the following dimensionless parameters:
 | (26) |
The bending modulus, 
(51) is of order 24 kBT (at 300 K), and we estimate a = 9 Å, J = 0.4 kBT a–1, 
= 0.7 kBT a–2, and g1 = 24 kBT a–2 (See near eqs. 6–9 for the estimates of a, J, and ). This leads to the following estimates:
 | (27) |
We assume g 1 kBT, which is of the order of the result from the Random Phase Approximation applied to a mixture of monomers (54). This leads to the estimates
 | (28) |
We control the quench depth by 
, where, within mean-field theory, 
leads to growth. Finally, the thickness lDOPC of a DOPC bilayer is 37 Å and that of SM is lSM = 51 Å, leading to 
, assuming ideal mixing l() = lDOPC + (1 – ) lSM.
Fig. 9 a shows the growth rate 
as a function of qa for a shallow quench depth 
and for the dimensionless measure of the hydrophobic surface tension 
from 0 to 25. For small , the membrane is stable. However, hydrophobic mismatch starts to drive phase separation upon increasing , with the dominant growing wavelength 
decreasing. Upon increasing the quench depth 
to –2 (Fig. 9 b), the membrane becomes unstable at = 0, and hence is below the critical point. There are two dominant growing wavelengths for this quench depth, one that increases and one that decreases upon increasing . The large wavelength instability is a perturbation of the ordinary spinodal instability after a quench in a two-phase mixture. The small wavelength instability is due to the coupling of concentration and thickness gradients; this coupling renders the theory unstable in the absence of higher order thickness and composition gradient terms, which stabilize the instability at a higher wavenumber. Upon further increasing the quench depth to 
(Fig. 9 c), there is only one dominant growing wavevector. The dominant growing wavenumber q*a versus for 
and increasing quench depths is shown in Fig. 9 d. The decrease in q*a (increase in *) upon increasing for equal to –2 and –3 is due to the fast thickness relaxation coupling to the composition, so that hydrophobicity perturbs the ordinary phase separation. The increase in q*a (decrease in *) upon further increasing for all quench depths shown is due to the competition between hydrophobicity and bending effects. The hydrophobic coupling induces a short wavelength instability (see Eq. 10), stabilized by the bending energy which is higher-order in wavenumber.
Next we will ensure that our mesoscopic and continuum models produce consistent results. Obviously these models address different time regimes and so a comparison between, for example, Figs. 1 and 9 is not possible. However, both models should accurately predict the existence of mixed and demixed membrane morphologies. To investigate this, we have determined the versus –f''0 phase diagram (Fig. 10) from the continuum model, where increasing –
corresponds to increasing the quench depth and therefore . As can be observed in Fig. 10, an increase in quench depth leads to a decrease in the value of required to induce phase separation. The same qualitative behavior is displayed by the lower phase boundary in the phase diagram determined numerically (Fig. 4). Note that the lower phase boundary in Fig. 4 corresponds to the case where finite size effects are unimportant. Obviously, this is the relevant phase boundary for comparison with Fig. 10 since finite size effects are also unimportant here. For quenches where 
phase coexistence occurs regardless of the hydrophobicity (Fig. 10). The discontinuity observed is due to the mean-field nature of our continuum model. We expect that fluctuations will smooth out the discontinuity because growing modes will become coupled.
Next we investigate the effect of decreasing the thickness relaxation time via increasing for and = 25 (Fig. 11 a) and 
and = 15 (Fig. 11 b). Changing , i.e., the system kinetics, does not effect the range of unstable modes or appreciably effect the dominant growing wavelength *. However, * does become more unstable upon increasing . Hence, a faster thickness response reduces the resistance to coarsening, allowing faster coarsening.
In summary, we have used a continuum model to investigate both the instability and the resultant initial domain growth of a membrane due to hydrophobic mismatch. We have found that the characteristic size for domain formation is a nonmonotonic function of the dimensionless measure of hydrophobic surface tension for the quench depths equal to –2 and –3. Firstly, domains increase in size upon increasing . Upon further increasing , domain sizes decrease for all quench depths investigated. Next, we investigated the effects of decreasing the thickness relaxation time via increasing . We found that upon an increase in , there was no appreciable change in the dominant growing wavevector. However, an increase in did allow faster coarsening.
|
CONCLUSIONS
|
|---|
TOP
ABSTRACT
INTRODUCTION
(b). Varying 
versus qa for different 
.
, 
and 
.