THE JOURNAL OF CHEMICAL PHYSICS 143, 243124 (2015)

Correlated lateral phase separations in stacks of lipid membranes
Takuma Hoshino,1,2,3,a) Shigeyuki Komura,1,3,b) and David Andelman2,3,c)

1

Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University,
Tokyo 192-0397, Japan2
Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University,
Ramat Aviv, Tel Aviv 69978, Israel
3
Kavli Institute for Theoretical Physics China, CAS, Beijing 100190, China(Received 14 August 2015; accepted 20 October 2015; published online 11 November 2015)
Motivated by the experimental study of Tayebi et al. [Nat. Mater. 11, 1074 (2012)] on phaseseparation of stacked multi-component lipid bilayers, we propose a model composed of stacked
two-dimensional Ising spins. We study both its static and dynamical features using Monte Carlosimulations with Kawasaki spin exchange dynamics that conserves the order parameter. We showthat at thermodynamical equilibrium, due to strong inter-layer correlations, the system forms a
continuous columnar structure for any finite interaction across adjacent layers. Furthermore, the phaseseparation shows a faster dynamics as the inter-layer interaction is increased. This temporal behavior
is mainly due to an effective deeper temperature quench because of the larger value of the criticaltemperature, Tc, for larger inter-layer interaction. When the temperature ratio, T/Tc, is kept fixed, thetemporal growth exponent does not increase and even slightly decreases as a function of the increased
inter-layer interaction. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4934984]

I. INTRODUCTION

Biological membranes are constructed out of two monolayers (leaflets) arranged in a back-to-back configuration.
They are mainly composed of phospholipids but contain also
other molecules such as cholesterol, glyco-sugars, and proteins.1 In living organisms, these membranes can form multilamellar stacks known as lamellar bodies.2 Examples of suchhighly folded membranous structures are thylakoid membranes of photosynthetic cyanobacteria or plant chloroplasts,
and stratum corneum of human skin. Since multilamellar structures can combine single membrane functions in series, they
offer possibilities for novel applications in photonics and as
bio-sensors.
Over the last decade, many studies have been performedon artificial giant unilamellar vesicles (GUVs) composed ofternary mixtures of saturated lipid such as sphingomyelin,
unsaturated lipid such as DOPC (1,2-dioleoyl-sn-glycero-3phosphocholine), and cholesterol.3,4 By decreasing temperature, these ternary mixtures undergo a lateral phase separation,
where a liquid-disordered (Ld) phase coexists with a liquidordered (Lo) one. It is known that the Lo phase is rich insaturated lipid and cholesterol, while the Ld phase is rich in
the unsaturated lipid.
In a recent experimental study, Tayebi et al.5 reported
that a stack (typically composed of several hundred layers) ofmulticomponent lipid bilayers with phase-separated domains
exhibits inter-layer columnar ordering. Using ternary mixturesof sphingomyelin, DOPC, and cholesterol, it was observedthat domains in stacked bilayers align one on top of the other,
thereby forming an uninterrupted columnar ordering acrossa)hoshino-takuma@ed.tmu.ac.jpb)komura@tmu.ac.jpc)andelman@post.tau.ac.il

0021-9606/2015/143(24)/243124/9/$30.00

hundreds of bilayer membranes. Such a cooperative multilayer epitaxy was attributed to the interplay between intralayer domain growth and inter-layer coupling. The formationof columnar structures in stacked bilayers is important because
it allows for electrical currents and transport processes to passthrough many transmembrane channels in a cooperative andefficient manner. Other possible applications of the columnar
ordering can be as templates for membrane protein crystallization, which is necessary for X-ray structural analysis ofmembrane proteins incorporated in bilayers.
As far as the dynamics of phase separation in stacks
of membranes is concerned, the temporal evolution of theaverage inplane domain size, R, was shown to obey a powerlaw growth, R ∼ t α with α ≈ 0.455.5 This exponent is larger
than the value obtained using GUVs with a single bilayer, forwhich the reported experimental value is α ≈ 0.28 ± 0.05.6
Hence, Tayebi et al. concluded that the inplane domain growthin each of the bilayers of the stack is faster, as compared to thedomain growth in GUVs.
In a subsequent paper,7 a model based on regular solution
theory, which takes into account the inter-lamellar coupling of
inplane phase-separated domains, was proposed. The calculated phase diagram was presented in terms of intra-layer andinter-layer coupling parameters and contains three different
regions: (i) a “one-phase” region in which the system does notexhibit phase separation; (ii) a “two-phase” region in which
two phases coexist and domains in different layers along thenormal z-direction are completely aligned and have the same
composition in the various layers, and (iii) a “multi-phase”
region in which there are unaligned inplane domains withdifferent compositions in the different layers. According toRef. 7, the transition line between the “two-phase” and “multiphase” regions strongly depends on the number of layers in thestack which was varied up to ten layers.

143, 243124-1© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-2

Hoshino, Komura, and Andelman

Being motivated by these works,5,7 we investigate thecorrelation between lateral phase separation in a stack of multilayer membranes using a spin model called the stacked twodimensional (2d) Ising model. This is the simplest model
to describe a stack of binary membranes composed of two
types of lipids. The model is the same as the anisotropic
three-dimensional (3d) Ising model for a finite stack in the zdirection. The important difference between the two models
is that in the former the order parameter (magnetization) ineach layer is conserved. This requirement is based on theexperimental fact that the A/B lipid composition in each layer
almost does not change during experimental times.
In our model, we study the thermodynamical equilibrium features using Monte Carlo (MC) simulations. The main
reason that we performed MC simulations rather than analyzing the mean-field free energy describing phase separation(as studied in Ref. 7) is to allow us to investigate the role ofthermal fluctuations on the inter-layer domain correlation instacked membranes. We show that the domains in each layer
are correlated along the vertical z-direction, for any finite valueof the inter-layer interaction is positive, i.e., J ′ > 0. Hence, thesystem is either in a one- or two-phase state in equilibrium,
and in our model the “multi-phase” state is not obtained in
the thermodynamic limit of infinite lateral size, as long as theinter-layer coupling J ′ > 0. As anticipated, it is found that thephase-transition temperature, Tc(J ′), increases as a function of
the inter-layer interaction parameter.
We also investigate the dynamics of phase separation at
fixed temperature T in the two-phase coexistence region. We
show that the accelerated temporal behavior of the phase separation for the stack is mainly driven by the increase of thetemperature quench, ∆T = Tc(J ′) − T, because Tc(J ′) becomeslarger for larger J ′. However, if the ratio T/Tc(J ′) is kept fixed,
the dynamics of the phase separation is actually slower forlarger values of the inter-layer coupling, J ′.
In Sec. II, we describe the stacked 2d Ising model andreview the MC simulation method. In Sec. III, we present theequilibrium properties of the model and discuss the condition for domain columnar ordering. Section IV describes the
dynamics of domain growth for different values of the interlayer interaction, and it is compared with a previous theoreticalwork.
II. MODEL AND SIMULATION TECHNIQUE

In our simulations, we use the stacked 2d Ising model,
shown in Fig. 1(a). We consider a stack of two-component
lipid bilayer membranes composed of an A/B lipid mixture,
although the experimental systems often consist of ternary
lipid/cholesterol mixtures. This simplification does not affectthe essential feature of the lateral phase separation. Anothersimplification is that we treat only symmetric bilayers where
the composition of the two leaflets is identical. Hence, eachlipid bilayer having a finite thickness can be mapped intoa 2d Ising model with conserved magnetization, expressing
the fact that no lipid is allowed to exchange across layers.
The 2d Ising layers are stacked in the z-direction, and theyinteract with their two nearest-neighboring layers, as depicted
in Fig. 1(b).

J. Chem. Phys. 143, 243124 (2015)

FIG. 1. (a) Schematic illustration of a stack of binary membranes, taken
here as a stack of three bilayers in the z-direction. Each bilayer is composedof two identical leaflets containing saturated lipids (A, black) and unsaturated lipids (B, white). Saturated and unsaturated lipids typically form Lo
and Ld phases, respectively. As the lipid molecules are not allowed to exchange between different bilayers, their composition in each bilayer is fixed.
(b) The stacked two-dimensional (2d) Ising model. Here the bilayer structureof each membrane is neglected. Lipids A and B correspond to spin up (black)and spin down (white), respectively. J is the coupling parameter between
nearest-neighbor spins in the same layer, while J ′ is the coupling parameter
between spins belonging to two nearest-neighboring layers.

The Hamiltonian of this stacked and coupled 2d Isingsystem can be written as
H = −JSi, ρ Si, ρ′ − J ′Si, ρ Si+1, ρi,⟨ρ, ρ ′⟩−

µi Si, ρ ,

i, ρ

(1)

i, ρ

where up/down values of the spin, Si, ρ = ±1, at ρ = (x, y) in
the ith layer correspond to a lattice site occupied by an A or
B lipid, respectively. The coupling between nearest-neighbor
spins in the x y-plane (denoted by ⟨ρ, ρ ′⟩) is J, while thecoupling with the nearest-neighbor spins across layers in the zdirection is J ′. The physical origin of the inter-layer interaction
J ′ is primarily attributed to direct van der Waals attractive
interactions acting between neighboring bilayers.8 Other nonspecific interactions, such as electrostatic and/or hydrationinteractions, can be taken into account through the second
virial coefficient and will affect the value of J ′ as well.9,10
Throughout this paper, we shall use the dimensionless ratio
defined by λ ≡ J ′/J as a measure of the inter-layer couplingstrength.
In the above Hamiltonian, µi is the external field (chemical potential), which fixes the average magnetization (A/B
composition) in the ith layer. Although µi can, in general, takedifferent values for different layers, we consider here the casewhere all of them are the same, µi = µ, fixing the same valueof lipid composition in all layers. This assumption holds alsoThis article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-3

Hoshino, Komura, and Andelman

for the dynamical states since we do not allow the lipids to beexchanged across different layers. The average order parameter(A/B composition) in the ith layer is denoted by S i, ρ , andthroughout this paper (except in Fig. 5(b)) we choose S i, ρ = 0,
which corresponds to a symmetric 1:1 A/B lipid mixture, i.e., at
the critical composition. This is equivalent to setting the valueof the chemical potential to zero, i.e., µ = 0.
The present model is related to the anisotropic 3d Isingmodel for a finite slab. The special case of λ = 1 corresponds
to the isotropic 3d Ising model, whereas for λ = 0 the stack
is composed of non-interacting 2d Ising layers. One interesting issue related to the anisotropic model, 0 < λ < 1, is the
crossover from 2d to 3d critical behavior11 that will be explored
below. We also note that the stacked 2d Ising model has beenstudied a great deal in connection with multilayer adsorption
phenomena on attractive substrates,12,13 but not in the context
of layers of binary mixtures with conserved magnetization(order parameter) as studied in this paper.
We investigate both the statics and dynamics of a stackof membranes based on the Hamiltonian presented in Eq. (1).
We employ MC simulations for classical Ising spins on a finite
L × L × L z lattice. Periodic boundary conditions are used in
all three directions. The spin configurations are updated using Kawasaki exchange dynamics14 in order to conserve themagnetization in each layer. This is based on the experimental
fact that lipids almost do not exchange across different layers.
Hence, their A/B inplane composition is fixed during experimental times.
The MC simulations presented here are performed in thefollowing way. At each MC trial step, a site on the 3d lattice and one of its nearest neighbors in the same layer arechosen at random. If the two spins are alike, a new site is

J. Chem. Phys. 143, 243124 (2015)

again chosen at random. This process is repeated until two
unlike nearest neighbor spins are found. Then, the probabilityof exchanging the two spins is determined by the standard
Metropolis algorithm.15 If the energy difference due to the spin
exchange becomes negative, i.e., ∆E < 0, we accept the exchange. Otherwise, we accept the exchange with a probabilityexp(−∆E/T), where T is the temperature and the Boltzmann
constant, k B, was set to unity.
In one Monte Carlo step (MCS), this procedure is repeated
L × L × L z times. The MC simulations are carried out by
annealing the temperature gradually from an initial infinite
temperature for which the spin configurations are completely
disordered and uncorrelated. The first 105 (or in some cases
up to 106) MCS is discarded in order to reach thermal equilibration. Furthermore, to avoid correlations between different
equilibrated configurations, measurements are taken every 20
MCS, and we averaged over 105 independent system configurations, in order to obtain sufficient statistics.
In order to investigate the phase separation dynamics, wemonitor the domain coarsening as a function of time (MCS) at
a constant temperature below Tc. An example of a typical time
evolution of phase separation is presented in Fig. 2 for λ = 0.1,
T/J = 1.63, and L = L z = 64, where six snap-shots are shown
from 102 MCS till 107 MCS. For clarity purposes, only theboundaries between domains of spin up (rich in lipid A) andspin down (rich in lipid B) are shown. In the initial time steps,
the phase separation occurs inplane, and the domains coarsen
without much out-of-plane coupling (due to the rather smallvalue of λ = 0.1). As time evolves, the inplane coarsening is
also followed by out-of-plane columnar ordering, where thelipid A (and lipid B) rich domains are highly correlated along
the z-direction. This is clearly seen for the fully equilibrated

FIG. 2. Time evolution of phase separated domains in the stacked 2d Ising model at different MC steps for λ = 0.1 and T /J = 1.63. The other parameters
are S i, ρ = 0 and L = L z = 64. For presentation purposes, only the interfaces of domain boundaries are shown, and the two different sides of the interfaces are
represented by green and brown. The system is fully equilibrated after about 107 MCS.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-4

Hoshino, Komura, and Andelman

J. Chem. Phys. 143, 243124 (2015)

configuration occurring after about 107 MCS (last snap-shot).
Here, the two color boundaries represent the two sides ofthe domain boundaries (while the inside of the domain is notshown). The boundaries look like extended interfaces separating inplane domains that are vertically connected along the zdirection, in agreement with experiment.5
III. STATIC PROPERTIES OF THE STACKED DOMAINS

In order to determine the phase-transition temperatureand obtain the corresponding phase diagram, we compute thespecific heat per lattice site defined as)
1 1 ( 2⟨H ⟩ − ⟨H⟩2 ,
(2)
c= 2
2
L Lz T
where H is given by Eq. (1) and ⟨· · · ⟩ indicates an ensemble
average. We note again that the above specific heat is calculated at constant magnetization (corresponding to constantlipid concentration in our model) of each layer. In our simulations, the ensemble average is taken by averaging over independent equilibrium spin configurations as explained in Sec. II.
For a given system size and dimensionless ratio λ, we calculatec as function of the dimensionless temperature T/J. Such adependence of c on T/J is presented in Fig. 3(a) for several
lateral system-sizes, L, and for λ = 0.1, L z = 8, recalling thatL z is the number of layers of the 3d stack.
For each system size, we associate the peak position of thespecific heat with the apparent critical temperature, Tc(L, λ),for a system with a finite size, L. Finite-size scaling analysis isthen performed in order to determine the critical temperature
for a slab of a finite L z layers in the thermodynamic limit(L → ∞). In Fig. 3(b), we plot Tc(L, λ = 0.1) as a function of
1/L for the same parameters as in (a). The plotted data are fittedwith the following finite-size scaling assumption:
Tc(L, λ) = Tc(λ) + aL −1/ν ,(3)

where Tc(λ) = Tc(L → ∞, λ) is the infinite system criticaltemperature for a given λ, a is a non-universal prefactor, andν is the 2d critical exponent for the correlation length in thex y-plane. We set ν = 1 in our analysis, following the work by
Pham Phu et al.,16 who performed extensive MC simulations
on magnetic Ising films (with λ = 1).17 We choose this 2d
critical exponent for the fitting because it was shown16 thatthe 2d character of the film is dominant even for L z = 13. The
extrapolated critical temperature for L → ∞ obtained from
Fig. 3(b) is Tc(λ = 0.1)/J = 2.85. We repeat this procedure fordifferent values of the inter-layer interaction parameter in therange of 0 ≤ λ ≤ 1 and determine the corresponding criticaltemperature, Tc(λ). We note that the value ν = 1 provides agood fitting for all the λ values examined.
Somewhat surprisingly, finite-size effects in the z-direction are much weaker as compared to those in the lateral
direction. This is shown in Fig. 4, where we plot c as a function of T/J when (a) λ = 0.1, (b) λ = 0.5, and (c) λ = 1
for a different number of layers, L z = 4, 8, 12, 16, while thelateral size L = 48 is kept fixed. For all λ values studied here
(0.1 ≤ λ ≤ 1), the observed peak position, T/J ≈ 2.65 in (a),3.30 in (b), and 4.10 in (c), is almost independent of L z , at least
for L z ≥ 8. This means that, in our model with a fixed imposed

FIG. 3. (a) Specific heat per lattice site, c, as function of the dimensionlesstemperature T /J , for different lateral system-sizes L = 16, 24, 32, 40, 48. The
other parameters are S i, ρ = 0, λ = 0.1, and L z = 8. For each system size,
the peak position of c is identified with an effective “phase transition” temperature. (b) Finite-size scaling analysis of the phase-transition temperature,
Tc(L)/J for λ = 0.1. The apparent phase-transition temperature is plotted as afunction of 1/L. The solid line is the fit given by Eq. (3) with ν = 1 (see text).
The extrapolated value for the critical temperature is Tc(λ = 0.1)/J = 2.85.

magnetization (A/B composition) in each layer, the correlationin the z-direction is very strong due to the cooperative behaviorof domains in different layers.
For fully equilibrated configurations, as shown in Fig. 2
after 107 MCS, the domains are highly connected verticallyalong the z-direction, from the bottom layer to the top one.
This is also shown in Fig. 5 in which the columnar structure ofdomains in different layers is clearly shown. Hence, the correlation length in this direction exceeds L z , and the constraint of
fixed magnetization (A/B composition) in each layer induces
a strong structural correlation in the z-direction even though
the inter-layer interaction J ′ is smaller than the intra-layerinteraction J (λ ≤ 1). A more quantitative argument for thedomain connectivity will be given later. Because the number ofThis article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-5

Hoshino, Komura, and Andelman

J. Chem. Phys. 143, 243124 (2015)

FIG. 4. Specific heat per lattice site, c, as function of the dimensionless temperature T /J , for different systems sizes L z = 4, 8, 12, 16 for (a) λ = 0.1, (b) λ = 0.5,
and (c) λ = 1. The other parameters are S i, ρ = 0 and L = 48. The observed peak position, T /J ≈ 2.65 in (a), 3.30 in (b), and 4.10 in (c), is almost independentof L z , at least for L z ≥ 8.

layers, L z , barely affects the MC results as shown in Fig. 4 for
λ = 0.1, 0.5, and 1, most of the simulations were done using
L z = 8, which is sufficiently large in our case to observe theasymptotic behavior of L z → ∞. For the anisotropic 3d Isingmodel without any constraint of conserved magnetization,
as previously studied in Ref. 11, a very weak system-size
dependence of the apparent critical temperature was observedby measuring the planar susceptibility.
The results of finite-size scaling analysis are shown in
Fig. 6, where we plot Tc as a function of λ. The criticaltemperature interpolates between the 2d and 3d Ising results,
toTc2d < Tc(λ) < Tc3d; the exact value in 2d (corresponding√
λ = 0) is known to be Tc2d/J = 2/ ln(1 + 2) ≈ 2.269 forsquare lattices,18 and the numerical estimate in 3d (corresponding to λ = 1) is Tc3d/J ≈ 4.511 for cubic lattices.19 These twolimits are recovered in our simulations and are seen in Fig. 6
for λ = 0 and 1, respectively. Although a more detailed
λ-dependent scaling behavior of Tc(λ) was previously discussed in the limit of very small λ,11,20 we shall generalize theargument for the anisotropic case of finite λ, 0 ≤ λ ≤ 1. When
T < Tc(λ), the stack undergoes a phase separation, and theinplane domains rich in lipid A (spin up) are interconnectedalong the z-direction, bridging between adjacent layers andforming large connected domains of the same average composition. The same feature also occurs for the B-rich domains.
Such a behavior can be clearly observed in Fig. 5.

In order to monitor quantitatively the degree of interconnectivity of domains in different layers, we define thefollowing quantity:
δ2 =


2
1 * 1 
+
S
−
S
,
i, ρ
i, ρ
L2 ρ , L z i-(4)where the average is taken over equilibrated MC configurations
as explained above. This quantity can be cast also as1 (Si, ρ − S i, ρ )(S j, ρ − S j, ρ ) ,(5)
δ2 = 2 2
L L z ρ i, j
and represents a special “magnetic susceptibility,” where thecorrelations are taken only along the z-direction and thenaveraged laterally in each of the planes. When the domains
are connected along the z-direction, the summation overdifferent i-layers will produce a large value of δ, while δis small if the domains are uncorrelated across the layers evenfor T < Tc(λ). In Fig. 7, we plot δ2 as a function of T/J
for different values of λ, while fixing L = 16 and L z = 8.
Notice that even for λ as small as 0.05 (blue diamonds), δ2
tends to increase as the temperature decreases below Tc(λ).
This means that the domains are connected in the z-direction
once the phase separation takes place. On the other hand,domains are independent and uncorrelated only when the interlayer interaction is extremely small, i.e., λ ≤ 0.001 in Fig. 7.

FIG. 5. Time evolution of phase-separated domains in a stacked 2d Ising model of eight layers, L z = 8, at different MC steps for (a) S i, ρ = 0 and (b) S i, ρ = 0.4.
The other parameters are λ = 0.1, T /J = 2.0, and L = 256.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-6

Hoshino, Komura, and Andelman

J. Chem. Phys. 143, 243124 (2015)

FIG. 8. Schematic representation of phase separated domains in a stack ofmembranes. Black and white domains are rich in A and B lipids, respectively. Two extreme cases are shown; (a) domains are fully connected in thez-direction, (b) domains are arranged at random and are disconnected.

semi-infinite domains: one composed by the A lipid (spin up)
and the other by the B lipid (spin down), as shown schematically in Fig. 8. When the domains are fully correlated in thez-direction, as in Fig. 8(a), the total free energy of the stackconsists of the contributions
FIG. 6. The phase-separation temperature, Tc(λ)/J , at the critical composition, as a function of the interaction parameter λ for symmetric A/B mixtures,
S i, ρ = 0. The system is in a phase-separated state below the solid line and ina one-phase state above the line.

The situation is found to be marginal when λ = 0.01 (redtriangles) because δ2 then slightly deviates from zero at low
temperatures.
Based on our MC results, we conclude that in the thermodynamic limit, L → ∞, domains will always be connected
for any finite inter-layer interaction, J ′ > 0. We give now asimple argument supporting this conclusion and show thatin the limit L → ∞ but with a finite number of layers, L z ,the domains in different layers are uncorrelated only whenJ ′ = 0 (λ = 0) is strictly obeyed. For the symmetric A/B case
(S i, ρ = 0), each layer will eventually phase separate into two

FIG. 7. The out-of-plane domain connectivity, δ 2, defined in Eq. (4), as afunction of the dimensionless temperature T /J , for different values of λ = 0,0.001, 0.01, 0.05, 0.1, 0.5, 1.0. The other parameters are S i, ρ = 0, L = 16,
and L z = 8. The transition temperatures for different λ values are indicatedby arrows. The value of δ 2 becomes larger when domains are correlated along
the z-direction between different layers. This increase in δ 2 is observed forlower temperatures and larger λ.

Fcon = −J ′L z L 2 + Fintra,(6)

where Fintra accounts for the intra-layer interactions. On theother hand, when the inplane domains are completely random
and disconnected, as sketched in Fig. 8(b), the total free energy
is dominated by an entropy contribution of arranging a random
stack of A and B domains along the z-direction,
Fdis = −T L z ln 2 + Fintra,(7)

with the same Fintra as before because this term is common forboth free energies. By comparing Eqs. (6) and (7), the thresholdinter-layer interaction, (J ′)∗, separating the two states, is givenby(J ′)∗ =T ln 2.
L2(8)

Notice that (J ′)∗ depends on L but not on L z . For finite temperatures, it vanishes in the thermodynamic limit of L → ∞.
Hence, this simple scaling argument suggests that domains arealways connected in the z-direction for any finite value of J ′.
Therefore, for all λ > 0, in the phase-separated region (below
the critical temperature) presented in Fig. 6, domains should
always form interconnected structures along the z-direction.
As shown in Eqs. (6) and (7), the internal energy scales with
L 2, while the entropy due to the random stacking of domains
does not depend on L. Hence, the entropic effect can neverovercome the internal energy in the thermodynamic limit andleads to the stability of the columnar structure. This conclusion
is not in agreement to that of Tayebi et al.,7 who claimed that
there is a “multi-phase” state in which domains are not aligned
and have different compositions even in thermodynamical
equilibrium.
In the simulations, (J ′)∗ can be finite due to finite-size effects. For instance, if the temperature is chosen to be
T/J = 1 in Fig. 7, the threshold value for L = 16 can be
estimated as λ ∗ = (J ′)∗/J ≈ 2.7 × 10−3. Since λ = 10−2 (redtriangles in Fig. 7) exceeds this threshold, the correspondingδ2 takes larger values at low temperatures. Moreover, the very
weak finite-size effects along the z-direction are consistent
with the lack of L z -dependence of (J ′)∗ in Eq. (8).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-7

Hoshino, Komura, and Andelman

J. Chem. Phys. 143, 243124 (2015)

IV. DYNAMICS OF PHASE SEPARATION

We address now the effects of inter-layer interaction on
the dynamics of phase separation as the system converges
towards its thermal equilibrium state. Under the assumptionthat scaling laws can be applied, the average domain size R
increases according to a temporal power-law: R(t) ∼ t α .4 For
2d systems for which the total domain area is conserved, theaverage domain size R is inversely proportional to the totalinterface length ℓ, i.e., R ∼ ℓ −1.21,22 This can easily be seen
because R and ℓ are related by ℓ = 2πnR ∼ nR, where n is
the number of domains, and the total area of all domains,
A = πnR2 ∼ nR2, is a conserved quantity. Hence, within thescaling hypothesis, the total interface length (in 2d) should
behave asℓ(t) ∼ t −α .

mention the scaling argument of Binder and Stauffer on phaseseparation dynamics of particles that undergo cluster reaction
and diffusion processes.24–27 Under the assumption that mostparticles that leave a cluster reimpinge on the same cluster atlater times, the diffusion coefficient D of a cluster of size R was
shown to scale as D ∼ R−(1+d), where d is the embedded spacedimension. If we further assume that the domain size R is the
only length scale in the system, the scaling relation for a simple
diffusion process is given by R2 ∼ Dt. This argument yields
the growth exponent to be α = 1/(3 + d). Hence, the predictedvalues from this scaling conjecture are α = 1/5 for d = 2 andα = 1/6 for d = 3.
Our simulation results, namely, α ≈ 0.24 for λ = 0 andα ≈ 0.14 for λ ≥ 0.2 compare favorably with this prediction.

(9)

In our stacked Ising model, we calculate the interface lengthin each of the layers and average it over different layers. Notethat the total interface length is proportional to the first term of
the Hamiltonian in Eq. (1), which enumerates the number ofbonds across the inplane A/B interface.
In Fig. 9(a), we plot the temporal evolution of the total interface length in 2d, ℓ(t) (and averaged along thez-direction), as a function of time measured in MC steps.
The temperature quench into the two-phase region is done
for a fixed temperature, T/J = 2.0 < Tc(λ), in order to mimic
the experiment that is conducted at fixed room temperature.
Several values of λ are studied, and the other parameters are
L = 256 and L z = 8, with averages taken over three independent MC runs. For each λ value, the scaling behavior of Eq. (9)is analyzed, and we extract the growth exponent α from thelate stage kinetics. We find that for λ = 0 (2d case), the growth
exponent has the smallest value of α ≈ 0.07, while for λ > 0,
it is a function of λ and increases up to α ≈ 0.28, as shown in
Fig. 9(b).
Although this result may explain the fact that the phaseseparation has an accelerated dynamics in stacked membranes
as compared to GUVs (isolated single membranes), we should
keep in mind that Tc(λ) increases as function of the inter-layer
coupling λ > 0, as shown in Fig. 6. As long as the final quench
temperature is fixed to T/J = 2.0, the temperature quench
depth defined by ∆T = Tc(λ) − T becomes larger as the valueof λ is increased. This may explain why the growth exponent
α becomes larger with increasing λ, for a fixed T-quench.
In order to have a better comparison between different
λ values, we evaluate in Fig. 10 the growth exponent in adifferent way. We now keep a constant quench ratio
T/Tc(λ) = 0.6, where T is the final quench temperature, andthe critical temperature Tc(λ) depends on λ, as shown in Fig. 6.
For these deeper temperature quenches (farther from Tc(λ)),
the estimated growth exponent is α ≈ 0.24 for λ = 0 (pure 2dcase), and 0.13 ≤ α ≤ 0.16 for 0.2 ≤ λ ≤ 1.0. Note that the
α-values are only weakly dependent on λ > 0.
Finally, we elaborate on the decreasing λ-dependence ofthe growth exponent α and show that this behavior is consistent
with the change in the dimensionality of the stack from 2d
to 3d. In general, the growth exponent associated with phaseseparation depends on the dimensionality.23 In this context, weFIG. 9. (a) The temporal evolution of the total interface length ℓ as a functionof time (MCS) for different values of λ = 0, 0.2, 0.4, 0.6, 0.8, 1.0, and fora temperature quench from the one-phase state (T → ∞) into the two-phasestate at T /J = 2.0. The A/B mixture is symmetric, S i, ρ = 0, L = 256, andL z = 8. The average over three independent MC runs is taken for each λvalue. The two dashed lines represent a power-law behavior with exponent
α = 0.07 and 0.28, which roughly bound the two limiting behaviors of the
λ-dependent exponent, α. (b) The domain growth exponent α as a functionof λ, as obtained from (a).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-8

Hoshino, Komura, and Andelman

FIG. 10. The temporal evolution of the total interface length ℓ as a functionof time (MCS) for different values of λ = 0, 0.2, 0.4, 0.6, 0.8, 1.0, and for atemperature quench from the one-phase state into the two-phase one, withfinal temperature satisfying T = 0.6Tc(λ). The A/B mixture is symmetric,
S i, ρ = 0, L = 256, and L z = 8. The average over three independent MC runs
is taken for each λ value. The two dashed lines represent a power-law
behavior with exponents α = 0.14 and α = 0.24. (b) The domain growth
exponent α as a function of λ, as obtained from (a).

The growth exponent decreases for finite λ because the systemcrosses over from 2d to 3d. This is due to the fact that thegrowing phase-separated domains are inter-connected along
the z-direction for λ > 0. It should be noted, however, thatthe absolute value of α obtained from the simulation is notuniversal but strongly depends on the quench depth as shownin Fig. 9. This explains why the above exponents are notin complete agreement with the simple scaling argument of
Binder and Stauffer.

J. Chem. Phys. 143, 243124 (2015)

stacked 2d Ising model given in Eq. (1). We use a Monte
Carlo simulation scheme with Kawasaki exchange dynamics
that conserves the order parameter in each layer, in order toinvestigate both equilibrium and dynamical features. Performing finite-size scaling analysis only in the lateral direction,
while keeping the stack thickness fixed (mimicking the experiment), we determine the phase-transition temperature, Tc(λ),
by changing the inter-layer interaction parameter λ = J ′/J. As
shown in Fig. 6, the phase-transition temperature interpolates
between that of the 2d and 3d Ising models.
One of our main conclusions is that domains in each oneof the layers are always interconnected along the z-direction,
forming a continuous columnar structure for any finite interlayer interaction J ′ > 0, as shown in Fig. 5. This domain
structure is in accord with the experimental findings for stacks
of few dozen to few hundred layers.5 However, the “multiphase” region in which there are unaligned inplane domains
with different compositions, as was predicted in Ref. 7, is
not found in our study at thermal equilibrium. On the other
hand, such a “multi-phase” state can be transiently observed
before the system reaches its fully equilibrated state, as can beobserved in Figs. 2 and 5.
We have also investigated the temporal evolution of
domain formation in the stacked 2d Ising model. Whenthe inter-layer interaction λ increases, the phase separation
appears to have an accelerated dynamics as can be seen by thelarger values of the growth exponent, α, shown in Fig. 9(b).
However, these larger α values are mainly due to an increasein the phase-transition temperature, Tc(λ), as function of λ;thus, a larger effective temperature quench, ∆T = Tc(λ) − T,for fixed T. When the final temperature quench T is fixed
relative to the phase-transition temperature as shown in Fig. 10
for T = 0.6Tc(λ), the growth exponent even decreases as the
λ value is increased. Our numerical findings for the growth
exponent α are different than the value of α ≈ 0.455, as
found in the experiment.5 One possible explanation for thisdiscrepancy can be the lack of hydrodynamic interactions in
our MC simulations.6
In this work, we have mainly discussed the case of
S i, ρ = 0, corresponding to the critical composition of the A/B
lipid mixture. Currently, we are investigating the dynamics
of phase separation for off-critical compositions, S i, ρ , 0[see Fig. 5(b)]. For such compositions, the phase-transition
temperature is smaller than the critical temperature. In thepresent simulations, the average A/B lipid composition (orderparameter of the Ising model) in each bilayer is restricted
to stay the same. In the future, we plan to study membranestacks where each layer has a different but fixed composition.28
Furthermore, since it is known from simulations that thepresence of a supporting solid substrate affects the dynamics ofmembrane domain growth,29 it will be of interest to incorporate
this substrate effect in future studies.
ACKNOWLEDGMENTS

V. CONCLUDING REMARKS

Motivated by recent works of Tayebi et al.,5,7 who studied
experimentally and theoretically the phase separation in stacks
of multi-component lipid bilayers, we have investigated theWe thank H. T. Diep, J.-B. Fournier, T. Kato, R. Okamoto,
H. Orland, M. Schick, P. Sens, S. Shimobayashi, and M.
Turner for useful discussions. T.H. acknowledges support ofTokyo Metropolitan University for an international student

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21

243124-9

Hoshino, Komura, and Andelman

exchange program and the hospitality of Tel Aviv University.
S.K. acknowledges support from the Grant-in-Aid for Scientific Research on Innovative Areas “Fluctuation and Structure”(Grant No. 25103010) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan, the Grant-inAid for Scientific Research (C) (Grant No. 24540439) fromthe Japan Society for the Promotion of Science (JSPS), and the
JSPS Core-to-Core Program “International Research Network
for Non-equilibrium Dynamics of Soft Matter.” D.A. acknowledges support from the Israel Science Foundation under Grant
No. 438/12 and the United States–Israel Binational Science
Foundation under Grant No. 2012/060.

J. Chem. Phys. 143, 243124 (2015)
9S. T. Milner and D. Roux, J. Phys. I (France) 2, 1741 (1992).
10S. Komura and D. Andelman, Europhys. Lett. 64, 844 (2003).
11K. W. Lee, J. Korean Phys. Soc. 40, L398 (2002).
12M. J. de Oliveira and R. B. Griffiths, Surf. Sci. 71, 687 (1978).
13R. Pandit, M. Schick, and M. Wortis, Phys. Rev. B 26, 5112 (1982).
14K. Kawasaki, in Phase Transitions and Critical Phenomena, edited by C.

Domb and M. S. Green (Academic Press, New York, 1972), Vol. 2.
15D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statis-

tical Physics (Cambridge University Press, New York, 2009).
16X. T. Pham Phu, V. Thanh Ngo, and H. T. Diep, Surf. Sci. 603, 109 (2009).
17T. W. Capehart and M. E. Fisher, Phys. Rev. B 13, 5021 (1976).
18L. Onsager, Phys. Rev. 65, 117 (1944).
19K. Binder and E. Luijten, Phys. Rep. 344, 179 (2001).
20L. L. Liu and H. E. Stanley, Phys. Rev. Lett. 29, 927 (1972).
21M. Laradji and P. B. S. Kumar, Phys. Rev. Lett. 93, 198105 (2004); J. Chem.

Phys. 123, 224902 (2005).
1B. Alberts, A. Johnson, P. Walter, J. Lewis, and M. Raff, Molecular Biology

22S. Ramachandran, S. Komura, and G. Gompper, EPL 89, 56001 (2010).

of the Cell (Garland Science, New York, 2008).

23A. J. Bray, Adv. Phys. 51, 481 (2002).

2G. Schmitz and G. Müller, J. Lipid Res. 32, 1539 (1991).

24K. Binder and D. Stauffer, Phys. Rev. Lett. 33, 1006 (1974).

3S. L. Veatch and S. L. Keller, Biochim. Biophys. Acta 1746, 172 (2005).

25A. B. Bortz, M. H. Kalos, J. L. Lebowitz, and M. A. Zendejas, Phys. Rev. B

4S. Komura and D. Andelman, Adv. Colloid Interface Sci. 208, 34 (2014).
5L. Tayebi, Y. Ma, D. Vashaee, G. Chen, S. K. Sinha, and A. N. Parikh, Nat.

Mater. 11, 1074 (2012).
6C. A. Stanich, A. R. Honerkamp-Smith, G. G. Putzel, C. S. Warth, A. K.

Lamprecht, P. Mandal, E. Mann, T.-A. D. Hua, and S. L. Keller, Biophys. J.
105, 172 (2013).
7L. Tayebi, A. N. Parikh, and D. Vashaee, Int. J. Mol. Sci. 14, 3824 (2013).
8J. N. Israelachivili, Intermolecular and Surface Forces (Academic Press,
New York, 2011).

10, 535 (1974).
26J. Marro, A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, Phys. Rev. B 12,2000 (1975).
27M. Rao, M. H. Kalos, J. L. Lebowitz, and J. Marro, Phys. Rev. B 13, 4328

(1976).
28K. Sornbundit, C. Modchang, N. Nuttavut, W. Ngamsaad, D. Triampo, andW. Triampo, J. Korean Phys. Soc. 63, 71 (2013).
29W. Ngamsaad, S. May, A. J. Wagner, and W. Triampo, Soft Matter 7, 2848

(2011).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
133.86.44.116 On: Fri, 13 Nov 2015 01:06:21