October 2017
EPL, 120 (2017) 18004
doi: 10.1209/0295-5075/120/18004

www.epljournal.org

Permeation through a lamellar stack of lipid mixtures
Takuma Hoshino1 , Shigeyuki Komura1 and David Andelman2
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
received 7 September 2017; accepted in final form 27 November 2017
published online 27 December 201787.16.D- – Membranes, bilayers, and vesicles
87.16.Uv – Active transport processes
PACS 05.40.-a – Fluctuation phenomena, random processes, noise, and Brownian motion
PACS
PACS

Abstract – We study material transport and permeation through a lamellar stack of multicomponent lipid membranes by performing Monte Carlo simulations of a stacked two-dimensional
Ising model in the presence of permeants. In the model, permeants are transported through thestack via in-plane lipid clusters, which are inter-connected in the vertical direction. These clusters
are formed transiently by concentration ﬂuctuations of the lipid mixture, and the permeation
process is aﬀected especially close to the critical temperature of the binary mixture. We showthat the permeation rate decays exponentially as a function of temperature and permeant lateralsize, whereas the dependency on the characteristic waiting time obeys a stretched exponentialfunction. The material transport through such lipid clusters can be signiﬁcantly aﬀected aroundphysiological temperatures.
c EPLA, 2017
Copyright Introduction. – Molecular transport and permeationin living systems are topics of growing interest in cosmetic and drug delivery applications [1]. Such processes,
characterized by selective transport of small ions andmacromolecules (biopolymers and proteins), often occur
through stacks of biomembranes. Examples of multilamellar structures within the cell are the Golgi apparatus
and mitochondria [2]. In thylakoid membranes of chloroplasts, the ordered array of one of the photosystems in astack of membranes is responsible for photosynthetic functions such as energy transfer and electron transport [3].
Moreover, it is known that different photosystems areheterogeneously distributed between the stacked and unstacked regions of thylakoid membranes [4].
On a much larger tissue scale, an important systemis the stratum corneum, which constitutes the outermost
layer of human skin [5]. It is known that the stratum corneum is composed of corneocytes and intercellularlipids that form lamellar structures [6]. Although stratum
corneum lipids form a gel phase with limited mobility at
physiological temperatures, the lipid tail-tail interface ofeach bilayer is in a liquid-like disordered state [7]. Cholesterol molecules incorporated in this liquid-like region candiffuse both translationally and rotationally, which allows

high overall cholesterol mobility. Hence such inter-leaflet
disordered regions can be regarded as heterogeneous fluid
sheets forming a multi-layered stack.
Using artificial stacks of multi-component lipid bilayers,
Tayebi et al. [8] reported that in-plane phase separationof lipid domains leads to an inter-layer columnar ordering between the domains. Such a strong vertical correlation between domains residing on adjacent membranes
can lead to a material transport, mediated through channel proteins that are preferentially incorporated into theselipid domains. However, it should be equally noted thata macroscopic phase separation usually does not occur inbiomembranes at physiological conditions, while thermal
fluctuations of local concentrations are always present insuch multi-component membranes [9–13]. The correlation
length characterizing concentration fluctuations increases
as one approaches the critical temperature of the lipid
mixture. For example, it was shown that sub-micrometer(about 50 nm) concentration fluctuations can take placein membranes at temperatures about 2–8 ◦ C above theircritical temperature [9,10].
In this letter, we investigate the permeation process
through a stack of two-component (saturated and unsaturated) lipid membranes that exhibit strong concentration

18004-p1

Takuma Hoshino et al.
fluctuations in the thermodynamically stable one-phaseregion (above the critical temperature). We regard themulti-layered binary membranes as a stack of twodimensional (2D) Ising model, and consider their interlayer correlations [14]. The permeable molecules residing
on one membrane can be transferred, within our model,
to the adjacent membrane through transiently connectedclusters bridging the two neighboring membranes. Using
Monte Carlo simulations, we investigate in detail the dependency of the permeation rate on temperature, permeant size, and the characteristic waiting time.
An important effect that is apparent close to the critical point is an enhanced permeability through stacks ofmembranes because the lifetime of transiently connectedclusters along the stack becomes longer due to criticalslowing-down. We demonstrate that the permeation rateincreases exponentially as the temperature approaches thecritical temperature from above, and conclude that thisrate is related to the in-plane concentration correlations.
Furthermore, the permeation rate is shown to decay exponentially as a function of the permeant lateral size.
We predict that concentration fluctuations in physiologicalconditions can play an important role for efficient materialtransport through multi-component membranes.
Model. – As shown in fig. 1, we consider a stack of lipid
membranes composed of a mixture of saturated (S) andunsaturated (U) lipids, modeled via the stacked 2D Isingmodel. Each lipid bilayer has a finite thickness and can bemapped into a 2D Ising model with conserved magnetization corresponding to the average S/U lipid composition.
The Hamiltonian of the membrane stack is [14]Si,ρ Si,ρ′ − J ′Si,ρ Si+1,ρ ,(1)
H = −Ji,ρ,ρ′ i,ρ

where the spin variable Si,ρ = ±1 is located at in-plane
position ρ = (x, y) of the i-th layer along the z-direction,and corresponds to a lattice site occupied either by an S or
U lipid, respectively. Furthermore, J is the in-plane coupling parameter between nearest-neighbor spins (lipids) in
the same layer, while J ′ is the coupling parameter between
spins (lipids) belonging to two nearest-neighboring layers,
originating primarily from direct van der Waals attractive
interactions [15]. Throughout this letter, we make use of
the inter-layer coupling strength, defined as the dimensionless ratio λ ≡ J ′ /J.
To simulate the behavior of permeants in a stack of lipid
membranes, Monte Carlo simulations are performed forthe equivalent Ising spin model, eq. (1). The simulations
are done on a finite 3D lattice of size L × L × Lz , with periodic boundary conditions in all three spatial directions.
The spin configurations are updated using Kawasaki exchange dynamics, in order to conserve the S/U lipid composition (the magnetization for the Ising model) in eachlayer [14]. Spin exchanges are only allowed for spin pairsthat belong to the same layer, because the time scale of lateral diffusion is much faster than that of out-of-plane lipid

Fig. 1: The stacked 2D Ising model used to model a stackof membranes when the inner structure of the membranes isignored [14]. Each membrane is composed of a mixture ofsaturated lipids (S) and unsaturated lipids (U). The S and U
lipids correspond to spin up (black) and spin down (white), respectively. The coupling parameter between nearest-neighbor
spins (lipids) in the same layer is J, and the coupling parameter between spins (lipids) belonging to two nearest-neighboringlayers is J ′ .

exchange between neighboring layers. The probability ofspin exchange is determined by the standard Metropolis
algorithm. The important difference between the present
stacked 2D Ising model and the ordinary anisotropic 3D
Ising model [16] is that in the former the magnetization(corresponding to the lipid composition) in each layer isconserved. The thermodynamical properties such as the
λ-dependent critical temperature of the stacked 2D Isingmodel were studied in detail in ref. [14].
Next, we address the permeation process in the model.
In addition to the spin variable that accounts for the S/Umixture, we introduce another variable to model the permeant molecules. We assume that each permeant has asquare shape of m = n2 lattice sites and lies flat within
one layer. In other words, its thickness is comparable tothat of the lipid bilayer and will not be further considered.
As schematically shown in fig. 2, we require the following conditions for a permeant to be transferred to
the neighboring layer below. i) Initially, all permeant
molecules are incorporated in a lipid cluster composed of
S or U lipids, located in the top layer, z = Lz . ii) The
Monte Carlo simulation is run continuously till the S (orU) lipid cluster in the upper layer lies above another cluster of the same S (or U) type in the adjacent layer below.
iii) The transient overlap between the two clusters along
the z-direction should persist over tw Monte Carlo steps(MCS). Once the conditions i)–iii) are satisfied, the permeant molecule is moved down to the adjacent layer below,
while preserving the permeant lateral position. iv) The
procedure is repeated, in a unidirectional way (no backward movement is allowed), till all permeants reach thebottom layer of the stack, z = 1.
The characteristic waiting time, tw in units of MCS,
introduced above is closely related to the diffusivity or
chemical affinity of the permeant (see “Discussion” below). If the cluster connectivity is destroyed within thewaiting time tw , either by lateral motion or disassembly

18004-p2

Permeation through a lamellar stack of lipid mixtures
16

12

8

4

0
0

1000

2000

3000

4000

-2

10

Fig. 2: (Color online) (a) Schematic illustration along thez-axis of a stack of lipid membranes in the presence of permeant molecules. Black and white lipids correspond to S andU lipids, respectively. (b) A schematic representation of thepermeation model (top view). A permeant molecule is allowed
to move from the z = i layer to the z = i − 1 layer below it,
if the following conditions are satisﬁed: i) A permeant is initially incorporated in a lipid cluster at the top layer, z = Lz .
ii) The lipid clusters should be connected in the z-direction
across neighboring layers. iii) If the formation of clusters in
the z-direction persists over waiting time tw (in units of MCS),
the permeant is moved from the z = i layer to the z = i − 1
layer below it. iv) The procedure is repeated, in a unidirectional way, till all permeants reach the bottom layer of thestack, z = 1.

of the lipid cluster, the permeant is forced to stay in thesame layer. Furthermore, the permeant lateral position isfixed because only the relative motion between permeants
and the lipid cluster is important. Note that we do nottake into account any interaction acting between permeant
molecules, which is justified in the dilute limit of permeants. We also assume that the above permeation process
is unidirectional (always permeates from top to bottom),because protein machines such as ion channels are responsible for directed transport of specific materials in biological membranes. In experiments, such transport processes
may be realized by applying external forces (e.g., electricfield or concentration gradient) to permeant molecules.
One MCS contains L × L × Lz spin updates, and thefirst 5000 MCS are discarded in order to reach thermal
equilibrium. After equilibration, we allow the permeants
to penetrate through the layers according to the algorithm

10

-3

3.0

3.2

3.4

3.6

3.8

4.0

Fig. 3: (a) Average z-position of permeants as a function oftime t measured in units of MCS after the initial 5000 MCS
are discarded for equilibration purposes. The parameter values
are τ = kB T /J = 3.0, m = 4, tw = 5 (in units of MCS) andλ = J ′ /J = 0.1. The system size is L × L × Lz = 48 × 48 × 16.
The data can be ﬁtted by a linear relation as given by eq. (2),with a permeation rate γ = 4.9 × 10−3 . The bars indicate thestandard deviations obtained by averaging over 20 independentruns. (b) The permeation rate γ as a function of the scaled
temperature τ for m = 4 and tw = 5. The exponential decay ofγ, as in eq. (3), is obtained with ﬁtting coeﬃcients, A = 0.9 and
B = −2.5. Note that the critical temperature value for λ = 0.1
is τc = 2.85 (shown by the arrow). This value is lower than
the temperatures used in the runs, which are always within
the one-phase region. The ﬁtting error is within the size of thesymbol.

introduced above. The same permeation procedure is repeated by running several Monte Carlo runs in order toimprove the statistics. The average S/U composition istaken to be at its critical value, i.e., 1:1 mixture of S
and U lipids. The inter-layer coupling strength is fixed
to λ = J ′ /J = 0.1, resulting in a critical temperaturekB Tc /J ≈ 2.85 [14], where kB is the Boltzmann constant. We emphasize that all Monte Carlo simulations
are conducted for temperatures larger than Tc , namely,
T > Tc (λ).
Results. – In fig. 3(a), we plot the average z-positionof permeants (measured in units of the layer number) as

18004-p3

Takuma Hoshino et al.
a function of time t (measured in units of MCS after
discarding the initial 5000 MCS) for scaled temperature,τ ≡ kB T /J = 3.0 > τc = 2.85, permeant size m = 4, andwaiting time tw = 5. The obtained result can be fitted bya linear relation between position z and time t,z = −γt + Lz ,

-3

(2)

where the slope γ corresponds to the permeation rate
—to be distinguished from the permeability coefficient introduced earlier in the literature [17,18]. Here the leastsquare fitting yields γ = 4.9 × 10−3 , meaning that it
takes about γ −1 ≈ 200 MCS for a permeant (of size
m = 4) to move down from the i-th layer to the adjacent(i − 1)-th layer. According to the proposed Monte Carloprocedure as described above, the permeant molecules are
moved down only to the adjacent layer below within asingle MCS. Hence the above permeation rate should bealways smaller than the inverse of the waiting time tw ,i.e., γ < 1/tw . The linear dependency in time shown in
eq. (2) is a consequence of the unidirectional permeation
process that is an important assumption within our model.
Moreover, γ is roughly a constant during the permeation
process, since the average time needed to move from thei-th layer to the (i − 1)-th one is the same for any i. On
the other hand, γ depends on the reduced temperatureτ , the permeant size m, and the waiting time tw , i.e.,γ = γ(τ, m, tw ). In the following, we shall examine thesedependencies in more detail.
We first focus on the temperature dependency, and plot
the permeation rate γ as a function of τ > τc on a semilog plot in fig. 3(b) for m = 4 and tw = 5. Note thatthe scaled critical temperature, τc = 2.85, is below therange of τ we examined. As can be seen in fig. 3(b), theobtained temperature dependency can be well fitted by an
exponential formγ(τ, m, tw ) = exp[−A(m, tw )τ + B(m, tw )],

-2

10

10

-4

10

3.0

3.2

3.4

3.6

3.8

4.0

3.0

3.2

3.4

3.6

3.8

4.0

-2

10

10

-3

Fig. 4: (Color online) (a) The permeation rate γ as a functionof the scaled temperature τ for diﬀerent permeant size, m = 1,
4 and 9. Other parameter values are tw = 5 and λ = 0.1. (b)The permeation rate γ as a function of the scaled temperature
τ for diﬀerent waiting times tw = 1, 3, 5 and 7. The otherparameter values are m = 4 and λ = 0.1. In both (a) and (b),γ decays exponentially with τ except when m = 1 in (a) and
tw = 1 in (b).

(3)

where the fitting coefficients A and B are functions of mand tw . Notice that A is mostly positive, indicating thatthe permeation rate γ decreases as the scaled temperature,τ > τc increases. In other words, the permeation rate increases exponentially as the temperature approaches thecritical temperature from above. This is an important result and holds quite generally in our simulation. Such atemperature dependency of γ can later be explained interms of the correlation length which characterizes thelipid cluster size.
The above results indicate that in stacked membranes,
lipid concentration fluctuations significantly affect the permeability, especially close to the critical point, τ  τc . As
τ → τc , the correlation length substantially increases andeven diverges. This leads to a large increase of the relaxation time of those clusters (critical slowing-down), which
can easily exceed the waiting time, tw . Since the clusters in adjacent layers are strongly correlated even for asmall value of the coupling parameter λ = J ′ /J = 0.1 [14],the permeation rate depends strongly on temperature asin eq. (3).
In fig. 4(a), we plot the permeation rate γ as a functionof the temperature τ for different permeant sizes m = 1, 4
and 9, while fixing tw = 5. Although γ hardly depends on
τ for m = 1, it exhibits strong temperature dependency
for m = 4 and 9. Notice that m = 1 corresponds to permeants that have the minimal size allowed in the present
simulation. On the other hand, for any fixed temperature
above τc , the decrease in γ is found to be significant as
the permeant size m is increased. This can be understoodsince it is more difficult for larger clusters to satisfy our
above-mentioned permeant transport conditions, i) andii). In fact, γ decreases nearly exponentially with m whenfixing the values of all other parameters. This result implies that the cross-section area of a permeant, m = n2 ,is the important factor in controlling the permeation because the area occupied by a permeant is kept constant
while it is transported to the neighboring layer. Such a

18004-p4

Permeation through a lamellar stack of lipid mixtures
To check the validity of the above-proposed scaling
form, we have plotted in fig. 5 the quantity ln γ as a function of −Aτ + B with A and B assumed to obey the scaling form as in eqs. (4) and (5) for m = 4, 9 and 16, and
tw = 2, 3, . . . , 9. A good data collapse can be seen forsmaller m values, while there is a systematic deviation
below the fitting line for m = 16 (triangles) as tw is increased. Note that even in the latter case of m = 16, a
linear relation between ln γ and −Aτ + B is maintained,
indicating that the scaling form of B in eq. (5) becomesinaccurate for larger m and tw .

Fig. 5: (Color online) Scaling plot of ln γ as a function of −Aτ +B for m = 4, 9 and 16, and tw = 2, 3, . . . , 9. See eqs. (3), (4),and (5). Diﬀerent symbols indicate diﬀerent m values (shownin the graph), while diﬀerent colors indicate diﬀerent tw values;
tw = 2 (black), 3 (red), 4 (light green), 5 (blue), 6 (magenta),
7 (deep green), 8 (cyan), and 9 (white).

size dependency is in accordance with the experimentallyobserved exponential dependency of the skin permeability
coefficient on molecular volume [18].
In fig. 4(b), we plot γ as a function of the temperature τfor different waiting times tw = 1, 3, 5 and 7, while fixing
m = 4. We remark that an exponential decay is found except for tw = 1, which is the smallest waiting time allowedin the model. Although γ is almost independent of τ for
tw = 1 and m = 4 due to these small values, γ generally
decreases with τ for larger values of m > 4 even for tw = 1(not shown here). Furthermore, γ is found to decrease as
tw increases for any fixed temperature. According to thepermeation condition iii), lipid clusters in two neighboring layers need to overlap vertically over tw time steps inorder to achieve the permeation. This implies that thepermeation rate γ should be a decreasing function of tw ,which is confirmed in our simulation. However, as discussed below, the dependency of γ on tw is not as simple
as its dependency on τ or m.
To discuss more quantitatively the behavior ofγ(τ, m, tw ) on m and tw , we analyze the functions A(m, tw )and B(m, tw ) of eq. (3). For the sake of simplicity, we assume that both A and B obey the following scaling forms:
A(m, tw ) ≈ a1 mα1 (tw )α2 + a2 ,(4)

B(m, tw ) ≈ b1 mβ1 (tw )β2 + b2 ,(5)with four exponents α1 , α2 , β1 , and β2 and four constants
a1 , a2 , b1 and b2 . The numerical fitting yields the following
values: α1 ≈ 1.03 ± 0.095, α2 ≈ 0.74 ± 0.067, and β1 ≈1.61 ± 0.16, β2 ≈ 1.37 ± 0.16. Then, the leading combined
dependency of γ on the three parameters, τ, m, and tw
can be written as γ ∼ exp[−m(tw )3/4 τ ], where only the
A-part dependency is written explicitly. This functional
dependence on the three parameters is the main result of
this letter.

Discussion. – The permeation rate γ, as is obtained from our simulation, decays exponentially withthe permeant size m and temperature τ . On the other
hand, the dependency on the waiting time tw is describedby a stretched exponential with an exponent α2 ≈ 3/4.
Our results are consistent with a stretched exponential relaxation associated with dynamic heterogeneity that canbe found in glassy systems. In the trapping model [19],a stretched exponential behavior with an exponent 3/5
was predicted for a dynamical correlation function in 3D.
Our obtained value α2 ≈ 3/4 is slightly larger, and thephysical meaning of the precise α2 value needs further
clarifications.
We point out that the temperature dependency of thepermeation rate in eq. (3) is analogous to that of the inplane correlation function of an Ising spin system. For
the one-phase region (T > Tc ), the correlation function is
given by
Si,ρ Si,ρ+r  ∼ exp(−r/ξ),(6)
where r = |r| and ξ is the correlation length. Close to
the critical point, the correlation length as a function oftemperature scales as ξ ∼ |τ − τc |−ν , where ν is the corresponding critical exponent. In our previous work [14], wehave shown that the critical behavior of the stacked Isingmodel can be described by a 2D Ising model even in thepresence of a coupling between adjacent layers. Hence, theexponent ν should take the exact 2D value of ν = 1 [20], as
has been also confirmed in our supplementary simulation(not shown here). Hence we have ξ ∼ |τ − τc |−1 .
Therefore, the in-plane correlation length is inversely proportional to the reduced temperature, i.e.,ξ ∼ 1/(τ − τc ), and the correlation function in eq. (6) decays exponentially with temperature. It is reasonable toexpect that the temperature dependency of the permeation rate is essentially the same as that of the correlation
function, because only the permeants incorporated in lipid
clusters can be transferred to the adjacent membrane.
Another finding of ref. [14] is that the correlations in
the z-direction are very strong because of the cooperative behavior of domains in different layers. This featurearises because of the constraint that the lipid compositionin each layer is strictly conserved. As a result, below thecritical temperature (τ < τc ), the system forms a continuous columnar structure for any finite interaction λ > 0across adjacent layers.

18004-p5

Takuma Hoshino et al.
Although so far we have discussed concentration fluctuations above the critical temperature (τ  τc ), we consider that the correlations in the z-direction are always
strong enough, especially close to the critical temperature.
Hence, we expect that the in-plane lateral correlation is agoverning factor for the temperature dependence of thepermeation rate γ. In the present work, we have studied
only the case of λ = 0.1. We think that this relatively
small coupling parameter is sufficient to have strong correlations in the z-direction. To understand the generalbehavior, however, it is necessary to investigate the permeation process for different λ-values.
As the lateral correlation length ξ increases close to
the critical temperature, the lifetime of lipid clusters also
increases due to critical slowing-down, and the typical relaxation time scales with ξ with a dynamical critical exponent [11]. The increased lifetime of lipid clusters at
criticality certainly enhances the permeation process, andshall be further investigated in future studies.
We remark that the obtained temperature dependencyof the permeation rate should be distinguished from thepreviously discussed Arrhenius-type behavior of the skinpermeability coefficient [17]. In experiments, the permeability coefficient was shown to obey the form Kp ∼exp(−Ea /kB T ), where Ea is an activation energy of solute
molecules or ions. Hence, Kp increases with temperature, and manifests an opposite trend when compared with
our result, eq. (3). In our work we have focused on theenhanced permeability due to the concentration fluctuations, which are very sensitive to the proximity to thecritical temperature Tc . If lateral heterogeneities in living biomembranes at physiological conditions correspondto critical fluctuations [9–11], we expect that the material transport through clusters should also be significantly
affected around physiological temperatures, as studied inthis letter.
A typical value of the permeability coefficient throughthe stratum corneum was measured [18] to be Kp ≈5 × 10−8 m/s. Assuming that the inter-membrane distance is about d ≈ 5 × 10−9 m, we can estimate a characteristic time scale for the permeation as d/K p ≈ 0.1 s,
which is fairly large. Although we cannot yet make a direct connection with the experiments, the waiting time
tw in our model should be comparable to this time scale.
From a microscopic point of view, the waiting time twis determined by the hydrophobic interactions between
the permeant molecule and the surrounding membrane
environment.
Finally, we note that the self-diffusion coefficient ofions through a single-component unilamellar vesicle wasshown to increase substantially near the lipid maintransition temperature [21]. In those studies, the phasetransition occurs between two states of lipid molecules:
an ordered gel phase (solid-like phase) and a disordered liquid crystalline phase (liquid-like phase). Furthermore, some models suggested [22–24] that the interfaces
between the gel and disordered-liquid domains of the lipid

molecules are responsible for the high molecular permeability. Our model differs from these models, because
we have focused on the enhanced permeability due toconcentration fluctuations in a multi-component lamellar
stack close to the critical temperature of the binary lipid
mixture.
Conclusions. – To summarize, we have performed
Monte Carlo simulations of a stack of binary lipid mixture
arranged in a multi-layered structure, and considered thetransport of permeant molecules through such a lamellar
stack. Within our model, permeants can be transported
to adjacent membranes through vertically connected lipid
clusters that are transiently formed by lipid concentration fluctuations in the one-phase above Tc . We havefound that the permeation rate decays exponentially withtemperature and with the permeant cross-sectional area,
whereas the dependency on the waiting time obeys astretched exponential behavior. Such an exponential dependence on temperature is analogous to the behavior of
the in-plane correlation function, and in accordance withthe experimentally observed dependency of the permeability coefficient on permeant size. Our results imply
that concentration fluctuations in physiological conditions
can play an important role for efficient material transport
through multi-component membranes.
It may be of interest to include in future works the excluded volume effect of the permeant molecules, as wellas the interaction between permeants, in order to broaden
the scope of the present model and connect it more directly
to permeation processes through multi-layered biologicalmembranes.
∗∗∗
We thank T. Kato, R. Okamoto and T. V. Sachin
Krishnan for helpful discussions. TH acknowledges
support by Grant-in-Aid for JSPS Fellows (Grant
No. 17J01643) from the Japan Society for the Promotion of Science (JSPS). SK acknowledges support by
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 (MEXT) of Japan, and by Grant-in-Aid
for Scientific Research (C) (Grant No. 15K05250) fromthe JSPS. DA acknowledges support from the Israel Science Foundation (ISF) under grant No. 438/12, the U.S.Israel Binational Science Foundation (BSF) under grantNo. 2012/060, and the ISF-NSFC joint research program
under grant No. 885/15. He also thanks Free University
Berlin, for its hospitality and the Alexander von Humboldt
Foundation for a Humboldt research award.
REFERENCES[1] Cevc G., Structure and Dynamics of Membranes,
Vol. 1A (Elsevier, Amsterdam) 1995.

18004-p6

Permeation through a lamellar stack of lipid mixtures[2] Alberts B. et al., Molecular Biology of the Cell (Garland
Science, New York) 2015.
[3] Iwai M., Yokono M., Kurokawa K., Ichihara A. andNakano A., Sci. Rep., 6 (2016) 29940.
[4] Borodich A., Rojdestvenski I. and Cottam M., Biophys. J., 85 (2003) 774.
[5] Hadgraft J. and Lane M. E., Phys. Chem. Chem.
Phys., 13 (2011) 5215.
[6] Hatta I., Nakazawa H., Obata Y., Ohta N., Inoue
K. and Yagi N., J. Phys.: Conf. Ser., 272 (2011) 012025.
[7] Das C., Noro M. G. and Olmsted P. D., Soft Matter,
10 (2014) 7346.
[8] Tayebi L., Ma Y., Vashaee D., Chen G., Sinha S. K.
and Parikh A. N., Nat. Mater., 11 (2012) 1074.
[9] Honerkamp-Smith A. R., Cicuta P., Collins M. D.,
Veatch S. L., den Nijs M., Schick M. and Keller
S. L., Biophys. J., 95 (2008) 236.
[10] Veatch S. L., Cicuta P., Sengupta P., HonerkampSmith A. R., Holowka D. and Baird B., ACS Chem.
Biol., 3 (2008) 287.
[11] Honerkamp-Smith A. R., Machta B. B. and Keller
S. L., Phys. Rev. Lett., 108 (2012) 265702.

[12] Ramachandran S., Komura S., Seki K. and Imai M.,
Soft Matter, 7 (2011) 1524.
[13] Komura S. and Andelman D., Adv. Colloid Interface
Sci., 208 (2014) 34.
[14] Hoshino T., Komura S. and Andelman D., J. Chem.
Phys., 143 (2015) 243124.
[15] Israelachivili J. N., Intermolecular and Surface Forces(Academic Press, New York) 2011.
[16] Lee K. W., J. Korean Phys. Soc., 40 (2002) L398.
[17] Scheuplein R. J. and Blank I. J., Physiol. Rev., 51(1971) 702.
[18] Potts R. O. and Guy R. H., Pharm. Res., 9 (1992)
663.
[19] Phillips J. C., Rep. Prog. Phys., 59 (1996) 1133.
[20] Onsager L., Phys. Rev., 65 (1944) 117.
[21] Papahadjopoulos D., Jacobson K., Nir S. and Isac
T., Biochim. Biophys. Acta, 311 (1973) 330.
[22] Nagle J. F. and Scott H. L., Biochim. Biophys. Acta,
513 (1978) 236.
[23] Cruzeiro-Hansson L. and Mouritsen O. G., Biochim.
Biophys. Acta, 944 (1988) 63.
[24] Yang L. and Kindt J. T., Langmuir, 31 (2015) 2187.

18004-p7