A Derivation of G_{\phi\phi}(r)

# Concentration fluctuations and phase transitions in coupled modulated bilayers

## Abstract

We consider the formation of finite-size domains in lipid bilayers consisting of saturated and hybrid lipids. First, we describe a monolayer model that includes a coupling between a compositional scalar field and a two-dimensional vectorial order-parameter. Such a coupling yields an effective two-dimensional microemulsion free-energy for the lipid monolayer, and its characteristic length of compositional modulations can be considered as the origin of finite-size domains in biological membranes. Next, we consider a coupled bilayer composed of two modulated monolayers, and discuss the static and dynamic properties of concentration fluctuations above the transition temperature. We also investigate the micro-phase separation below the transition temperature, and compare the micro-phase separated structures with statics and dynamics of concentration fluctuations above the transition.

## I Introduction

Biomembranes are two-dimensional (2D) fluids that separate the inner and outer environment of organelles in biological cells. Naturally occurring biomembranes consist typically of numerous lipid species, sterols, sugars and membrane proteins. According to the “lipid raft” hypothesis (1), some of the lipid components and/or proteins are incorporated into finite-size domains, which play an important role on cellular functions such as signal transduction processes. Recent experiments suggest that lipid rafts are nothing but dynamical molecular assemblies of 20 nm in size with finite lifetime in the order of 10–20 ms (2).

Being motivated by the raft hypothesis, a large number of investigations (3); (4); (5); (6); (7); (8); (9); (10); (11); (12); (13); (14) have been conducted to reveal the properties of artificial membranes consisting of lipid mixtures and cholesterol. Below the miscibility transition temperature, formation of micron-size domains were observed using fluorescent microscopy (3); (4); (5); (6). In some cases, rather than a macroscopic phase separation, domains with distinct size in the micrometer range have been reported (7); (8); (9). For example, various types of modulated (stripe or hexagonal) patterns have been found for multicomponent lipid and cholesterol mixtures (10). Above the miscibility transition temperature, even multicomponent membranes do not phase separate, and their concentration fluctuations around the homogeneous state can be investigated (11); (12); (13), in particular, close to the critical point, . Furthermore, it is interesting to note that critical concentration fluctuations have been observed in membranes extracted from living cells (14).

One of the main reasons that initiated the notion of lipid rafts in biomembranes is the existence of finite-size domains, rather than domains resulting from a macroscopic phase separation (15). Assuming that membranes are in equilibrium, the same question can be phrased in terms of whether the biomembrane state is above or below the miscibility transition temperature. In the high temperature one-phase region, the only relevant length scale is the correlation length associated with concentration fluctuations, and it diverges at the critical point,  (16). Below this temperature, there should be a physical mechanism suppressing domain coarsening in order to explain the existence of finite-size domains in equilibrium.

Yet another characteristic feature of biomembranes is that the lipid composition of the two leaflets (monolayers) constituting the bilayer is not the same (17). Moreover, such asymmetric monolayers are not independent but are coupled to one another. This was confirmed experimentally by investigating the phase separation of bilayers with different monolayer lipid compositions (18), or seen in simulations (19). One mechanism that leads to the coupling between the two leaflets is the lipid chain interdigitation occurring at the mid-plane of the bilayer (20); (21); (22), which may affect the domain size in asymmetric bilayers.

One possibility to account for such finite-size domains in lipid mixtures is to consider the special role played by “hybrid lipids” (23); (24); (25); (26). These lipids have one saturated hydrocarbon chain and another unsaturated one. They tend to be localized at boundaries of 2D domains and act as a line-active agent. Hence, a ternary lipid mixture consisting of saturated lipids with two saturated chains, unsaturated lipids with two unsaturated chains, and hybrid lipids can be regarded as a “2D microemulsion”, because it is analogous to the microemulsion phase consisting of oil, water and surfactant (27); (13). For microemulsions it is known that there is another length scale, in addition to the correlation length mentioned above, related to the size of water and oil micro-domains (28).

Using the microemulsion analogy, we suggest that an additional length scale can be responsible for the finite-size structures in biomembranes. We propose a model for bilayers consisting of two coupled monolayers that are both in a 2D microemulsion state, manifesting modulated phases. First, we describe a monolayer model that includes a coupling between the lipid composition and a 2D vectorial order-parameter. The model accounts for the line active nature of hybrid lipids. Next, we consider a bilayer in which two such monolayers are coupled through an inter-leaflet interaction. For the coupled bilayer, we discuss the static and dynamic properties of concentration fluctuations above , and investigate their micro-phase separation below  (29). Intermediate structures arise when two competing structures have different characteristic length scales. One of our important conclusions is that the micro-phase separated structures below reflect the static and dynamic properties of concentration fluctuations above .

This paper is organized as follows. In Sec. II, after describing how the microemulsion state is obtained for monolayers, a model for coupled modulated bilayers is presented. In Sec. III, we show the results of static and dynamic structure factors for the coupled bilayers above . In Sec. IV, we describe some results for competing micro-phase separation in coupled bilayers below . Finally, in Sec. V, further discussion and final remarks are presented.

## Ii Model

### ii.1 Modulated lipid monolayers

We consider a 2D microemulsion formation in a monolayer consisting of two types of lipids: a saturated lipid (denoted by “S”) and a hybrid one (“H”). As shown schematically in Fig. 1(a), the saturated lipid has two saturated chains, whereas the hybrid lipid has one saturated chain and an unsaturated one. It is well known (3); (4); (5); (6) that the resulting liquid-ordered () and liquid-disordered () phases are rich in saturated and hybrid lipids, respectively, as is depicted in Fig. 1(c). In the experiments mentioned previously, cholesterol is usually added as a third component, and is known to affect the area per headgroup of lipids (30); (31). However, since cholesterol has a strong preference for the phase and affects mainly the saturated lipid, we do not consider cholesterol explicitly in our model, and neglect its presence hereafter.

The local area fraction of saturated and hybrid lipids are defined as and , respectively, where is a 2D vector. Under the incompressibility condition, , and the only relevant order-parameter is the difference between the two lipid compositions, . The phenomenological free-energy that describes the 2D phase-separation between S and H lipids is given by a Ginzburg-Landau expansion in terms of the order-parameter :

 Fs[ϕ]=∫d2r[σ2(∇ϕ)2+τ2ϕ2+14ϕ4−μϕ]. (1)

Here is related to the line tension between monolayer domains, is the reduced temperature with respect to the critical point , and is the chemical potential that regulates the average value in the monolayer. Without loss of generality, the coefficient of the quartic term is set to be a positive constant by an appropriate rescaling of the position variable, .

Next we discuss the role of the hybrid lipid and its effect on the phase separation. We define a 2D lateral vector pointing from the unsaturated tail of the hybrid lipid toward its saturated one, as depicted in Fig. 1(a). We then introduce a coarse-grained 2D vectorial order-parameter field , as shown in Fig. 1(b). This vector is the spatial average of the vectors over areas large as compared with molecular size, but still small enough as compared with macroscopic scales. The phenomenological free-energy associated with the vectorial field (and due to the hybrid lipids) can be written as an expansion up to quadratic order in  (32):

 Fh[m]=∫d2r[K2(∇⋅m)2+a2m2]. (2)

The coefficient is the 2D elastic constant, while is taken to be positive so that is the stable homogeneous state. An additional term is allowed by symmetry, but is not considered at present for simplicity sake.

At the interface, hybrid lipids orient their saturated and unsaturated chains toward the and phases, respectively, thereby reducing the chain mismatch (25); (26) as shown in Fig. 1(c). In fact, the line activity of hybrid lipid was predicted (25); (26) to be more pronounced for binary mixtures of saturated and hybrid lipids as compared to ternary mixtures. Within our phenomenological approach, the role of hybrid lipids can be represented by a coupling between the lateral variation of and the vectorial field . To lowest orders, the coupling term in the free energy is given by

 Fc[ϕ,m]=−Γ∫d2rm⋅(∇ϕ), (3)

where is a positive coupling constant because the vector tends to orient towards the domains. The total monolayer free-energy is given here by the sum of the three terms introduced in Eqs. (1)–(3): . Notice that our model is valid in the weak segregation limit (close to ), because slow spatial variation of the order parameters is intrinsically assumed. We further note that a similar free energy functional using a vectorial order-parameter was proposed for 3D microemulsions (33).

The total free energy can be conveniently expressed using the 2D Fourier transform of

 Extra open brace or missing close brace (4)

where is a 2D wavevector, and similarly is the Fourier transform of . Minimizing with respect to , we obtain its optimum value as

 m(q)=iΓqa+Kq2ϕ(q), (5)

where . By substituting back Eq. (5) into the monolayer free energy , the minimized free energy is expressed as

 Fm[ϕ]= ∫d2q(2π)212[τ+σq2−Γ2q2a+Kq2]ϕ(q)ϕ(−q) +∫d2r[14ϕ4−μϕ], (6)

where the Fourier transform has been used only for the second-order terms in . Expanding the effective binary interaction for small , we obtain

 Fm[ϕ]≈ ∫d2q(2π)2[2Bq4−2Aq2+τ2]ϕ(q)ϕ(−q) +∫d2r[14ϕ4−μϕ], (7)

where two new parameters are defined

 B≡KΓ24a2,     A≡14(Γ2a−σ). (8)

When the coupling constant is small enough, i.e., (or ), the minimum of the Gaussian term (first integral) in Eq. (7) occurs at , and is a signature of a macroscopic phase separation (for ). On the other hand, when is large enough, i.e., (or ), the minimum occurs at a non-zero wavenumber given by

 q∗=√A2B=√a(1−aσ/Γ2)2K, (9)

indicating a potential micro-phase separation with modulation.

Taking the inverse Fourier transform, Eq. (7) can be expressed in position space as

 Fm[ϕ]=∫d2r [2B(∇2ϕ)2−2A(∇ϕ)2 +τ2ϕ2+14ϕ4−μϕ]. (10)

This is the so-called “2D microemulsion” free-energy for a monolayer composed of a binary lipid mixture. When , the negative gradient-squared term favors spatial modulations, while the positive Laplacian squared term with suppresses modulations. As mentioned above, Yamamoto et al. (25); (26) showed that the effective line tension between domains becomes negative for membranes consisting of saturated and hybrid lipids.

More generally, models based on equations similar to Eq. (10) have been used successfully in the past to describe modulated phases (34) arising in a variety of different biophysical and chemical systems such as Langmuir films (35), lipid membranes (36); (37) and diblock copolymers (38).

### ii.2 Coupled lipid bilayers

We consider two coupled modulated monolayers forming a bilayer as shown in Fig. 2. Each monolayer is a binary mixture of saturated and hybrid lipids. We define two local order-parameters for the two monolayers: and , depicted in Fig. 2. The coarse-grained free-energy functional is then written as (29):

 Fb[ϕ,ψ]=∫d2r [2B(∇2ϕ)2−2A(∇ϕ)2 +τϕ2ϕ2+14ϕ4−μϕϕ +2D(∇2ψ)2−2C(∇ψ)2 +τψ2ψ2+14ψ4−μψψ−Λϕψ]. (11)

The first five terms depend only on and its derivatives and describe the upper monolayer in Fig. 2 and its possible modulations. These are the same terms as in Eq. (10) for the single monolayer case. Similarly, the latter five terms describe the lower monolayer, where and are the two reduced temperatures, while and are the corresponding chemical potentials. The last term, , represents the coupling between the two leaflets with a coupling constant ,

We would like to comment on the physical origin of the proposed coupling term, , in Eq. (11). Note that this quadratic term is invariant under exchange of . When , this term can be obtained from a term such as , which represents a local energy penalty when the upper and lower monolayers have different compositions (20); (21). For mixed lipid bilayers, such a coupling may result from the conformational confinement of the lipid chains, and hence would have an entropic origin (20). By estimating the degree of lipid chain interdigitation, the order of magnitude of the coupling parameter can be evaluated (22). In general, can also be negative depending on the specific coupling mechanism (39). However, as the sign of does not affect our result, it is sufficient to consider the case. The other possible higher-order coupling terms which are allowed by symmetry are , and . However these terms do not affect the properties of concentration fluctuations in any essential way, and will not be considered hereafter.

## Iii Concentration fluctuations above Tc

Using the bilayer free-energy, Eq. (11), we obtain the static and dynamic structure factors, which describe the properties of concentration fluctuations for two coupled monolayers in their respective one-phase region (above ). We shall closely follow the formulation of Ref. (40) in which the coupled macro-phase separation was discussed for bilayers.

### iii.1 Static structure factor

The spatially varying and can be written as and , respectively, where and are the spatially averaged monolayer compositions, and and describe the deviations from their average values. In thermal equilibrium, and satisfy

 τϕϕ0+ϕ30−μϕ−Λψ0=0, τψψ0+ψ30−μψ−Λϕ0=0. (12)

Expanding the free energy, Eq. (11), and retaining the quadratic order terms in and , we obtain the Gaussian free-energy

 FG[δϕ,δψ]=∫d2r[2B(∇2δϕ)2−2A(∇δϕ)2+ϵϕ2(δϕ)2 +2D(∇2δψ)2−2C(∇δψ)2+ϵψ2(δψ)2−Λ(δϕ)(δψ)], (13)

where the notations and as well as Eq. (12) are used.

The static partial structure factor is , and because of the radial symmetry, , where . It is given by

 Sϕϕ(q)=2gψ(q)4gϕ(q)gψ(q)−Λ2, (14)

and similarly for and

 Sψψ(q)=2gϕ(q)4gϕ(q)gψ(q)−Λ2, (15)
 Sϕψ(q)=Sψϕ(q)=Λ4gϕ(q)gψ(q)−Λ2, (16)

and

 gϕ(q)=2Bq4−2Aq2+ϵϕ2, gψ(q)=2Dq4−2Cq2+ϵψ2. (17)

Since the structure factors diverge at , we see that the coupling parameter effectively shifts the critical temperature to lower values.

### iii.2 Decoupled leaflets (Λ=0)

When the two leaflets are decoupled, , it is sufficient to present results for only one of the two monolayers, say the one. From Eq. (14), the decoupled structure factor is simply given by

 Sϕϕ(q)=12gϕ(q). (18)

This function has a peak at where has a minimum, indicating that monolayer fluctuations have a characteristic length scale describing their spatial modulations. The peak height is hence given by

 Sϕϕ(q∗ϕ)=1ϵϕ−A2/B, (19)

and diverges at . In addition, the structure factor in Eq. (18) decays as for large wavenumbers.

The correlation function , where , is obtained by the inverse 2D Fourier transform of Eq. (18) (see Appendix A for the derivation)

 Gϕϕ(r)=ξϕλϕ32πBRe[H(1)0(2πrλϕ+irξϕ)], (20)

where is the zeroth order Hankel function of the first kind, and “Re” stands for the real part. This correlation function contains two length scales; the first is the modulation periodicity

 λϕ2π=(Bϵϕ)1/42√ 1−γϕ, (21)

and the second is the correlation length

 ξϕ=(Bϵϕ)1/42√1+γϕ, (22)

where .

In Fig. 3, we plot Eq. (20) for certain parameter values. The correlation function has exponentially decaying oscillations as a function of the distance , similar to 3D microemulsions (28). When , both and are finite, and the corresponding phase is called the structured-disordered phase. The peak of the structure factor occurs at for , whereas for . The line at is the Lifshitz line. The modulation periodicity diverges for , and the monolayer transforms into a disordered phase for . The line of is called the disorder line and is not a phase transition line. On the other hand, the correlation length diverges for (hence, this is the critical point), and the ordered phase appears for  (28).

As shown in Fig. 4(a), the phase diagram of a decoupled bilayer is easily obtained by combining the phase sequences in terms of the two independent parameters and . The various phases are expressed by the binary combination of ordered (O), structured-disordered (S) and disordered (D) phases for each of the leaflets. The two disorder lines are shown in the figure by dashed lines, while the Lifshitz lines occur at and (not shown).

### iii.3 Coupled leaflets (Λ≠0)

When the two monolayers are coupled (), the corresponding phase diagram can be obtained by analyzing the poles in the complex plane of Eqs. (14)–(16). One example is shown in Fig. 4(b). The boundary between the OO and SS phases (solid line) is determined by the condition that all the structure factors diverge. In the DD phase, all the poles are pure imaginary, whereas in the SS phase the poles are complex, and at least one real pole exists in the OO phase.

In Fig. 4(b) we see that the asymmetric phases such as DO and SD are suppressed as compared with Fig. 4(a). This is because one of the leaflets induces modulations in the other leaflet due to the inter-leaflet coupling, . In addition, the OO and SS phases are expanded as compared to Fig. 4(a). As mentioned before, increasing the coupling effectively lowers the bilayer temperature. The Lifshitz lines (dotted lines) are obtained by plotting the main peak positions of Eqs. (14) and (15). Due to the coupling effect, they are tilted as compared to the case where the lines occur at and .

In Fig. 5(a) and (b), we plot the structure factors of the decoupled and coupled bilayers, respectively. As an illustration of the coupling effect with , we consider two different characteristic wavelengths . The heights of the two peaks are set to be equal by requiring and (see Eq. (19)). The peak height of at is increased (as in Fig. 5(b)) due to the coupling effect, whereas that of at is almost unchanged compared with the decoupled case. We also plot given by Eq. (16) which represents the cross correlation of fluctuations between the two monolayers. This quantity is proportional to the coupling constant and its peak position is essentially determined by that of at .

The peak position of denoted as for the coupled case () is obtained numerically and plotted in Fig. 6(a) as a function of and , which are, respectively, the peak positions of and for the decoupled case. We find that for the coupled case, the value of is almost equal to that of the smaller of . In Fig. 6(b), the peak height of at is plotted as a function of the coupling parameter and the ratio of the two characteristic wavenumbers for the specific case of , as the parameters are fixed to and . The peak height increases as is increased because the temperature of the coupled bilayer is effectively lowered. To see this clearly, we have plotted in Fig. 6(b) the phase-transition line (dashed line) separating the structured-disordered and ordered phases. The peak height of also increases as approaches unity, which is the case where and completely overlap.

### iii.4 Dynamic structure factors

The dynamical fluctuations in composition, and , depend on time and are now considered for coupled modulated bilayers. Since both and are conserved quantities for each monolayer, the time evolution equations are given by

 ∂δϕ(r,t)∂t=Lϕ∇2δFGδ(δϕ)+νϕ(r,t), ∂δψ(r,t)∂t=Lψ∇2δFGδ(δψ)+νψ(r,t), (23)

where and are the kinetic coefficients taken to be constants, and represent Gaussian white noise, satisfying and

 ⟨νi(r,t)νj(r′,t′)⟩=−δijLi∇2δ(r−r′)δ(t−t′), (24)

where and now indicates the average over space and time.

The Fourier transform of fluctuations in both space and time is

 δϕ(q,ω)=∫d2rdtδϕ(r,t)e−i(q⋅r−ωt), (25)

and a similar Fourier transform is used for and . The dynamic structure factors such as are given by

 Sij(q,ω)=αij(q)ω2+[ω+(q)]2+βij(q)ω2+[ω−(q)]2, (26)

where the explicit expressions of and are given in Table 1 and the characteristic frequencies are:

 [ω±(q)]2=12[ω2ϕ(q)+ω2ψ(q) ∓√[ω2ϕ(q)−ω2ψ(q)]2+4LϕLψq4Λ2ω2ϕψ(q)], (27) ω2ϕ(q)=4L2ϕq4[gϕ(q)]2+LϕLψq4Λ2, (28) ω2ψ(q)=4L2ψq4[gψ(q)]2+LϕLψq4Λ2, (29) ωϕψ(q)=2q2[Lϕgϕ(q)+Lψgψ(q)], (30)

and and are defined in Eq. (17). In the more special case of , reduce to a simpler form

 ω±(q)=Lq2[ gϕ(q)+gψ(q) ∓√(gϕ(q)−gψ(q))2+Λ2]. (31)

The intermediate structure factors depend explicitly on time and are obtained by taking the inverse Fourier transform in :

 Sij(q,t)=∫dω2πSij(q,ω)e−iωt, (32)

and from Eq. (26), we obtain

 Sij(q,t)=αij(q)2ω+(q)e−ω+(q)t+βij(q)2ω−(q)e−ω−(q)t. (33)

Hence the decay of concentration fluctuations is described by a sum of two exponentials with two characteristic times . However, when the two monolayers are decoupled (), the and structure factors decay with a single exponential characterized by a decay rate and (see Eqs. (28) and (29)), respectively.

For the decoupled case (), the decay times and are plotted as a function of in Fig. 7(a), with the same parameters as those in Fig. 5. The plots show a shoulder reflecting the characteristic structure at wavenumbers and . Notice that the larger the initial length scale, the longer the decay time. For the coupled case with , we plot in Fig. 7(b). Due to the coupling effect, the two decay times split into a larger and smaller one, . The larger one, , exhibits two shoulders, while the smaller one, , has a shoulder between the two characteristic wavenumbers. The coupling affects the decay time of the structure corresponding to the smaller wavenumber (larger length), similar to the effect seen in Fig. 5 for the static structure factor.

## Iv Coupled modulated phases with different periodicities q∗ϕ≠q∗ψ

In our previous paper (29), we investigated the phase behavior of coupled modulated bilayers by using the free energy, Eq. (11), below . When the two monolayers have the same preferred periodicity, , we obtained the mean-field phase diagram exhibiting various combinations of modulated structures such as stripe (S) and hexagonal (H) phases. In some cases, the periodic structure in one of the monolayers induces a similar modulation in the second monolayer. Moreover, the region of the induced modulated phase expands as the coupling parameter becomes larger.

However, when the preferred periodicities in the two leaflets are different, , it is difficult to obtain the phase diagram because the free energy densities cannot be obtained analytically. We then have to rely on numerical simulations (29); (41) to solve the time evolution of the two coupled order-parameters as explained below.

For the dynamics of coupled modulated bilayers below , we use the following time evolution equations (29)

 ∂ϕ(r,t)∂t=Lϕ∇2δFbδϕ, ∂ψ(r,t)∂t=Lψ∇2δFbδψ, (34)

where the bilayer free-energy, , is given by Eq. (11). Below the effect of thermal fluctuations is less important and the noise terms have been omitted in the above equations. Hereafter, the kinetic coefficients and are set to unity for simplicity.

We solve numerically the above 2D equations using periodic boundary conditions. Each run starts from a homogeneous state with a small random noise around the average compositions and . Time is measured in discrete time steps, and corresponds to a well equilibrated system. In all our simulations, the reduced temperature parameters are fixed to be . The characteristic wavenumber in the monolayer is fixed as by setting , while the periodicity in the monolayer is varied and the condition is used as before. In the following, we shall consider only the two coupled stripe phases for . However, the combination of stripe and hexagonal phases leads to a rich variety of complex patterns (41).

When the characteristic wavenumbers of the two decoupled monolayers are different, , the two coupled modulated structures cannot apparently match each other. The frustration between the two different periodicities is due to the inter-leaflet coupling and affects their morphologies. In Fig. 8, we show one example of the time evolution of a coupled micro-phase separation for ( and ) and . The spatial patterns of , , , and the Fourier transformed pattern, , are presented for time steps of and . Starting from the isotropic state, the monolayer forms first stripes (a). Then, as the monolayer starts to segregate at larger , it simultaneously chops the stripes into smaller sections (b), and the monolayer transforms into a finger-like patterns (c). Reconnection of the stripes takes place after a long time of annealing, and a pattern of alternating fingers in the monolayer is finally obtained (d). From the time evolution of the 2D Fourier patterns of the monolayer, it is apparent that the intermediate structures are characterized by the two length scales of ratio .

In Fig. 9, we show the spatially modulated patterns at of the two coupled and monolayers with different periodicities ( and as before) for different values of the coupling parameter and . Notice that Fig. 9(b) is the same as Fig. 8(d). In the weak-coupling case (), the two monolayers exhibit two independent stripe morphologies with characteristic periodicities (called the “independent” morphology). Here the two stripes essentially do not affect each other. As the coupling constant is increased (), stripes with a finger-like structure appear in the monolayer, while the stripe morphology in the monolayer is almost unaffected (called the “intermediate” morphology). In the Fourier transformed pattern of , we clearly see that the structures with two different characteristic wavenumbers are coexisting. This result is also in accord with the properties of the static structure factors, , having two different characteristic lengths, as shown in Fig. 5. For a larger coupling parameter (), very similar patterns are obtained for and , and almost coincide with one another (called the “coincident” morphology). It should be noted that the structure with the larger wavelength dominates when the coupling is large enough, and the sequence of morphological changes shown in Fig. 9, as increases, is rather typical.

In order to quantify the three morphologies (independent, intermediate, coincident), we calculate the spatial average of the product of the two compositions

 ⟨ϕψ⟩=1A∫d2rϕ(r,t)ψ(r,t), (35)

where is the total system area. In Fig. 10, we plot at (sufficient for equilibration) for various combinations of and (for the case ). The morphology of the obtained patterns is marked on the figure by circles, crosses and squares. The values of are small for “independent” (circles) structures, while they become larger for the “intermediate” (crosses) and “coincident” (squares) morphologies. Although the morphological changes are gradual and do not represent a sharp transition, the intermediate patterns appear roughly for . The region of intermediate structure expands as increases, and the patterns coincide for . We note that although the morphology cannot be solely determined by the quantity , the behavior of is similar to that of the peak values of the cross correlation, , presented in Fig. 6(b).

For the quantitative argument concerning the micro-phase separation dynamics in each monolayer, we have also calculated the two self quantities:

 ⟨ϕ2⟩=1A∫d2rϕ2(r,t), (36)

and

 ⟨ψ2⟩=1A∫d2rψ2(r,t). (37)

In Fig. 11, we plot the square root of these quantities as a function of , with the same parameters as those used in Fig. 9. In all studied cases, the modulation of the monolayer having a larger wavenumber (), grows faster than that of the monolayer. We also remark that the structure formation in the monolayer is accelerated for larger , whereas that of the monolayer is almost unchanged. According to the linear stability analysis of Eq. (34), the initial growth rates of the unstable modes are essentially equivalent to the decay rates of the concentration fluctuations given in Eq. (31). This is consistent with Fig. 7, where the characteristic growth time () for larger is smaller than that for smaller . The growth rates increases with (see Fig. 11(a)) because the coupling effectively reduces the temperature and enhances the phase transition. The faster the decay of concentration fluctuations, the faster the structure formation.

## V Discussion and final remarks

In this paper we present a model for coupled modulated lipid bilayers. We start by considering a monolayer consisting of a mixture of saturated and hybrid lipids, and propose a phenomenological model that includes a coupling between the lipid composition and a 2D vectorial field. This coupling arises from the line-active nature of the hybrid lipid, which adjusts its tail orientation in order to reduce the line tension. Minimization of the monolayer free-energy with respect to the vectorial field yields a 2D microemulsion exhibiting modulated phases. The characteristic wavelength of modulation is determined by the monolayer coefficients and , Eq. (8), reflecting molecular properties of lipid mixture. We then construct a model for lipid bilayers comprised of two modulated monolayers that influence each other through an inter-leaflet coupling.

Based on the model, we study concentration fluctuations of bilayers above , and calculate their static structure factors. The calculated phase diagram for coupled bilayers shows that the extent of ordered and structured-disordered phases become larger as compared to the decoupled case. When the two monolayers with different preferred wavenumbers are coupled, (say ), the peak height of occurring at smaller -numbers becomes larger as compared to the decoupled case, whereas the peak height of occurring at larger -numbers almost does not change. Namely, the inter-leaflet coupling strongly affects the compositional modulation in each monolayer. Furthermore, the inter-leaflet coupling has a clear signature on the cross structure factor, , as well as on the dynamics of concentration fluctuations. By calculating the intermediate structure factor, , we show that concentration fluctuations generally exhibit a double exponential decay with two decay rates, . One of the decay times () exhibits two shoulders at wavenumbers describing the monolayer compositional modulations.

For membranes below , we studied the micro-phase separation of a coupled modulated bilayer. When the two monolayers have different modulations, , we obtained numerically a variety of complex patterns. The initial growth rates of the unstable modes are identical to the decay rates of the concentration fluctuations.

As mentioned in Sec. I, the special character of concentration fluctuations in our model may explain the finite-size domains (“rafts”) in biological membranes. For ordinary binary mixtures above , the only length scale of the disorder phase is determined by the correlation length, and close to , this length becomes large. But within our 2D microemulsion model, there is another length scale characterizing the modulations as given by Eq. (21). This second length scale may also explain the finite-size domains observed in some experiments as a result of micro-phase separation in the low temperature phase.

A related model based on a microemulsion picture was recently proposed by Schick (42), who considered a coupling between curvature and compositional asymmetry between the two leaflets, resulting in a 2D microemulsion. Although Schick’s model as well as ours share the microemulsion viewpoint, the origin of the physical mechanism is different. In Schick’s model the coupling between composition asymmetry and curvature gives rise to domains with different spontaneous curvature. Our model, on the other hand, assumes that for flat monolayers, the line-active nature of hybrid lipids is solely responsible for the microemulsion formation (see Eq. (3)). Hence, our proposed physical mechanism for the bilayer coupling, as discussed in Sec. II.2, is different. One of the consequences of our coupling is that the domains residing on the two leaflets are in register, unlike the prediction of Ref. (42).

The present work is concerned with the analysis of equilibrium properties and relaxation dynamics towards the equilibrium state. The existence of finite-size domains may also be explained by non-equilibrium lipid transport between the cell interior and the membrane. Such a mechanism was considered (43); (44); (45); (46) through a coupling between the membrane and an outer reservoir of lipid or cholesterol. Similarly to our model, these works explain the appearance of the finite-size domains as a result of micro-phase separation.

Asymmetry between lipid composition in the inner and outer leaflets of biomembranes has a very deep significance and is closely related to the cell biological functions. For instance, the breakdown of such compositional asymmetry is related to programmed cell death (apoptosis) (47). For living cells, this asymmetry is maintained by an enzyme called “flippase” that actively transports the lipids between the two leaflets (48). The half-time of lipid composition due to flip-flop motion is measured using time resolved small angle neutron scattering (SANS) (49), and is estimated to be several hours at physiological conditions.

Although membranal signal transduction is important for various biological functions, its dynamical properties are still not so well understood. As our theory offers an explanation for the size and dynamics of lipid domains, we hope that it and similar models will contribute in the future towards the understanding of functional processes in biomembranes.

###### Acknowledgements.
We thank T. Kato, S. L. Keller, S. Ramachandran, M. Schick and K. Yamada for useful discussions. YH acknowledges a Research Fellowship for Young Scientists No. 215097 from the JSPS. This work was supported by Grant-in-Aid for Scientific Research (grant No. 21540420) from the MEXT of Japan. SK acknowledges support by the JSPS Core-to-Core Program “International research network for non-equilibrium dynamics of soft matter”. DA acknowledges support from the Israel Science Foundation (ISF) under grant No. 231/08 and the US-Israel Binational Foundation (BSF) under grant No. 2006/055.

## Appendix A Derivation of Gϕϕ(r)

In this Appendix we present the derivation of of Eq. (20). The real-space correlation function is given by the 2D inverse Fourier transform of Eq. (18):

 Gϕϕ(r) =∫d2q(2π)2Sϕϕ(q)eiq⋅r =14πB∫∞0dqqJ0(qr)q4−(A/B)q2+ϵϕ/4B, (38)

where is the zeroth-order Bessel function. We use the relation where and are the zeroth-order Hankel functions of the first and second kind, respectively. Then, the integral in Eq. (38) is written as where

 Ii=∫∞0dqqH(i)0(qr)q4−(A/B)q2+ϵϕ/4B, (39)

with .

In order to evaluate the above integral, we performed the integration in the complex plane by replacing with the complex variable . The integrand has four poles at

 zj = (ϵϕ/B)1/42(±√1−γϕ±i√1+γϕ) (40) = ±2πλϕ ±i1ξϕ,

with and and have been defined in Eqs. (21) and (22). These poles are located in the quadrants 1, 2, 3 and 4, off the real and imaginary axes. For the integral , we integrate along a quarter-circle contour in the first quadrant in an anti-clockwise direction. The contour radius is taken to infinity (Fig. 12), and using the residue theorem we obtain

 ∫∞0dxxH(1)0(xr)x4−(A/B)x2+ϵϕ/4B +∫0∞dyiyH(1)0(iyr)y4+(A/B)y2+ϵϕ/4B=2πiRes|z1, (41)

where “Res” stands for the residue. Similarly, for the integral , we integrate along the contour of the quarter-circle of infinite radius in the fourth quadrant in the clockwise direction:

 ∫∞0dxxH(2)0(xr)x4−(A/B)x2+ϵϕ/4B +∫0−∞dyiyH(2)0(iyr)y4+(A/B)y2+ϵϕ/4B=−2πiRes|z4. (42)

Combining Eqs. (41) and (42), and further using the relation , we obtain

 I=π2√ϵϕ/B√