Enthalpic and Entropic Contributions to Interleaflet Coupling Drive
Domain Registration and Anti-registration in Biological Membrane
Akshara Sharma
∗
Department of Physics, Indian Institute of Science-Bangalore,
C. V. Raman Road, Bangalore, Karnataka 560012, India
Aniruddha Seal
School of Chemical Sciences,
National Institute of Science Education and Research,
Bhubaneswar, Khurda, Odisha 752050, India
Sahithya S. Iyer
†
and Anand Srivastava
‡
Molecular Biophysics Unit, Indian Institute of Science-Bangalore,
C. V. Raman Road, Bangalore, Karnataka 560012, India
(Dated: September 29, 2021)
Biological membrane is a complex self-assembly of lipids, sterols and proteins organized as a
fluid bilayer of two closely stacked lipid leaflets. Differential molecular interactions among its di-
verse constituents give rise to heterogeneities in the membrane lateral organization. Under certain
conditions, heterogeneities in the two leaflets can be spatially synchronised and exist as registered
domains across the bilayer. Several contrasting theories behind mechanisms that induce registra-
tion of nanoscale domains have been suggested[1–3]. Following a recent study[4] showing the effect
of position of lipid tail unsaturation on domain registration behavior, we decided to develop an
analytical theory to elucidate the driving forces that create and maintain domain registry across
leaflets. Towards this, we formulated a Hamiltonian for a stacked lattice system where site variables
encapsulate the lipid molecular properties including the position of unsaturation and various other
interactions that could drive phase separation and interleaflet coupling. We solve the Hamiltonian
using Monte Carlo simulations and create a complete phase diagram that reports the presence or
absence of registered domains as a function of various Hamiltonian parameters. We find that the
interleaflet coupling should be described as a competing enthalpic contribution due to interaction
of lipid tail termini, primarily due to saturated-saturated interactions, and an interleaflet entropic
contribution from overlap of unsaturated tail termini. We find that higher position of unsaturation
provides weaker interleaflet coupling. We also find points in our parameter space that allow ther-
modynamically stable nanodomains in our bilayer model, which we have verified by carrying out
extended Monte Carlo simulations. These persistent non-coalescing registered nanodomains close
to the lower end of the accepted nanodomain size range also point towards a possible “nanoscale”
emulsion description of lateral heterogeneities in biological membrane leaflets.
I. INTRODUCTION
The biological cell membrane contains a host of lipids
and other biomolecules such as proteins and glycans,
which interact dynamically to facilitate biological pro-
cesses like signal transduction and membrane protein
oligomerization. Lipids in model ternary and quater-
nary membrane systems can segregate into liquid ordered
(L
o
) and liquid disordered (L
d
) phases. Experimentally
observed phase separated domains from in vitro studies
range from micrometers to nanometers in their size[5–
12]. Apart from studies on model in vitro systems, seg-
regation of lipids into L
o
and L
d
phases has also been
∗
Also at Molecular Biophysics Unit, Indian Institute of Science-
Bangalore.
†
Currently at Department of Chemistry, University of Illinois
Chicago, 845 W. Taylor St., 4500 SES (MC 111), Chicago, IL
60607, U.S.A
‡
anand@iisc.ac.in
observed in in vivo membranes[13–16]. The functional
importance of this phase separation is exemplified dur-
ing signal transduction across cells, where the relay of
signal from the outer to inner leaflet and vice-versa is
highly dependent on the communication between the two
leaflets. Experimental studies have shown the impor-
tance of lipids in synchronising the reception and trans-
mission of messages across the bilayer[17–19]. T-cell and
B-cell receptor mediated immune response is a classical
example of systems where co-localization of outer and
inner leaflet ordered domains, i.e ”registration” of or-
dered lipids, plays an important role in lipid mediated
signal transduction[20–22]. Understanding the origin and
molecular driving forces giving rise to ”registration” or
”anti-registration” of the ordered domains in the two
leaflets can provide important insights into the related
physiologically critical processes.
Existing studies on the topic contain a variety of the-
ories, some of them conflicting with each other, on the
nature of the driving forces for domain registration. Ex-
change of cholesterol between ordered domains in the
.CC-BY 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted September 30, 2021. ; https://doi.org/10.1101/2021.09.29.462263doi: bioRxiv preprint
2
outer leaflet and the inner leaflet in asymmetric mem-
branes has been hypothesised to form and maintain reg-
istered ordered domains in the inner leaflet in some the-
oretical models[23–25]. A mismatch energy, contributed
when ordered and disordered domains overlap between
leaflets, has also been studied as a possible driving force
for domain registration[1]. This mismatch energy can be
understood in terms of the energy penalty required when
L
o
and L
d
domains interact between leaflets. In con-
trast, however, several studies claim membrane undula-
tions and domain bending stiffness to cause domain reg-
istration without any direct interleaflet coupling[3, 26].
There is also work done suggesting domain boundary
line tension as a driving force for registration and this
model also does not require explicit coupling of the lipid
bilayer leaflets[27]. These studies argue that the differ-
ence in splay rigidities between the L
o
and L
d
domains in
both leaflets causes preferential distribution of stiffer do-
mains in regions of the membrane with lower curvature
fluctuations. This, along with the energy gain from a
decrease of line tension at the domain boundaries, is sug-
gested as a possible explanation towards a driving force
for registration of both large and small domains. There
have been theoretical and computational studies on an-
other candidate for a possible driving force, hydropho-
bic mismatch, arising due to a difference in tail lengths
leading to differential membrane thickness between L
o
and L
d
domains[2, 16, 28–31]. Both domain formation
kinetics and registration dynamics have been studied us-
ing hydrophobic mismatch as a driving force, and it has
been theoretically shown to produce phase coexistence
metastable states that capture both registered and anti-
registered possibilities upon phase separation.
While there is agreement that there isn’t one single
driving force for this phenomenon, the complexity of the
problem comes from the number of hypotheses avail-
able that seem to explain domain registration or anti-
registration. Further additions have been made to this
list of candidates recently, of particular interest being a
computational study by S. Zhang and X. Lin[4] showing a
significant effect of the position of unsaturation along the
lipid tail on domain registration tendencies in otherwise
identical systems. The authors carried out MARTINI[32]
coarse-grained molecular dynamics (CGMD) simulations
of two systems that were a mixture of DPPC, Cholesterol
and an unsaturated lipid with different positions of un-
saturation, which they called D23 and D34. Fig.1 shows
the chemical structures of the lipids.
They observed that the system where the unsaturated
lipid had lower position of unsaturation showed registra-
tion and the other system showed anti-registration. They
also suggested that the interleaflet coupling was through
an attractive interaction in the interleaflet region between
the tail termini of the lipids.
We found this to be a very intriguing face of the prob-
lem, and decided to investigate the physical cause of
the position of unsaturation affecting domain registra-
tion characteristics, which we believed could provide im-
portant insights into the overall registration mechanism.
Besides testing the deductions made in the original work
that credited the enthalpic interaction in the interleaflet
region for driving registration, we also wanted to explore
the role of the competing entropic factors in the inter-
leaflet region. Our primary hypothesis is as follows: A
higher position of unsaturation would lead to reduced
configurational entropy of the lipid tails in the core of the
membrane leading to a reduction in interleaflet interac-
tion. On the other hand, a lower position of unsaturation
forces the tail terminus to explore the interleaflet region
(instead of bending towards the polar headgroups), which
would lead to not only a better enthalpic but also an en-
hanced entropic contribution to the interleaflet coupling.
The entropic aspect of interleaflet coupling is an impor-
tant and often overlooked aspect that we bring to light
in this work.
To test this, we formulate a Hamiltonian with tunable
parameters that captures the hypothesised enthalpic and
entropic contributions as a function of position of un-
saturation along the length of the lipid tail for a lattice
system representing the bilayer as a stack of two square
lattices. We used this Hamiltonian and conducted a pa-
rameter study using Monte Carlo simulations for each
point in parameter space for all the systems under inves-
tigation. Using these simulations, we try to capture the
most dominant factor affecting the registration character-
istics of the systems, as well as the effect of the change
in position of unsaturation for a system described by our
Hamiltonian. The schematic of our workflow is provided
in sec.I of the Supporting Information (S.I) as fig:S1.
We describe our lattice model of the membrane bilayer
in the ”Model” section below and provide simulation de-
tails, specifications of the studied systems, as well as de-
tails about our analysis tools in the ”Methods” section.
Following that, we report our observations and findings
in the ”Results” section and we put all our results in the
context of existing literature and provide a detailed per-
spective on the implications of our observations in the
”Discussion” section. We finally summarize our work
FIG. 1: The chemical structures of the DPPC, D23
and D34 lipids. Note that the D23 and D34 lipids differ
only in their position of unsaturation and do not occur
naturally.
.CC-BY 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted September 30, 2021. ; https://doi.org/10.1101/2021.09.29.462263doi: bioRxiv preprint
3
in the ”Conclusions” section with some perspectives on
position of unsaturation in lipids as a driving force for
domain registration.
II. MODEL
Following our hypothesis, we wrote a Hamiltonian for a
system with two coupled membrane leaflets that captured
entropic behaviour of the lipid tails as well as the posi-
tions of unsaturations in the lipid tails. For this, we de-
cided to model the system as a pair of stacked square lat-
tices, each square representing a membrane leaflet. Each
leaflet would be populated with lipids in a symmetric
population distribution. Our site variables for the lattice
Hamiltonian are carefully curated to capture the molec-
ular details of the lipid, which we describe below. Our
Hamiltonian can be broken into two independent parts,
one for the lateral interactions within a leaflet, and the
other for the interleaflet interactions.
H = H
lat
+ H
inter
(1)
Both H
lat
and H
inter
consist of an enthalpic and an en-
tropic part. H
lateral
is formulated as shown below where
the first term can be denoted as H
enth
lat
and the second
term as H
entr
lat
.
H
lat
=
1
2
X
<i,j>
X
k=0,1
X
α,α
0
αα
0
δ
α,ik
δ
α
0
,jk
+
V
S
k
B
T
4
X
k=0,1
X
i
n
α,ik
ln(φ
α,ik
)
+ n
α
0
,ik
ln(φ
α
0
,ik
)
(2)
Similarily, H
inter
is formulated as shown below, where
the first term can be denoted as H
enth
inter
and the second
term as H
entr
inter
.
H
inter
=
X
i
X
α,α
0
0
αα
0
δ
α,i0
δ
α
0
,i1
+ V
0
S
k
B
T
X
k=0,1
X
i
1 +
P
p
α
l
α
n
0
α,ik
ln(φ
0
α,ik
)
+
1 +
P
p
α
0
l
α
0
n
0
α
0
,ik
ln(φ
0
α
0
,ik
)
(3)
In the above eqns.2 and 3, the summation < i, j > is
a sum over orthogonal neighbours within a leaflet, and
the sum over i is a sum over the sites of a leaflet. The
δ
α,ik
are kronecker deltas that become 1 when the lipid
occupying site i in leaflet k is of species α. We used a
simple Ising-like enthalpic interaction between lipids in a
leaflet for the lateral part at site i as follows:
αα
0
= −V
αα
0
s
α,i0
s
α
0
,j0
, (4)
H
enth
lat
=
1
2
X
<i,j>
X
k=0,1
X
α,α
0
αα
0
δ
α,ik
δ
α
0
,jk
(5)
and for lipids opposite each other in the two leaflets
for the interleaflet part as shown:
0
αα
0
= −V
0
αα
0
s
α,i0
s
α
0
,i1
(6)
H
enth
inter
=
X
i
X
α,α
0
0
αα
0
δ
α,i0
δ
α
0
,i1
(7)
where V
αα
0
and V
0
αα
0
are the corresponding interaction
strength parameters that can be tuned and serve as a
design variable. The site variables s
α,ik
for species α at
site i in leaflet k = 0 or 1 is given by:
s
α,ik
=
l
α
X
c=1
S
α,c
l
α
(8)
which is the average S
CD
of the lipid species occupying
the site i in leaflet k, calculated from lipid trajectories
generated from all-atom molecular dynamics (AAMD)
simulations described ahead. As can be seen, the forms
in eqns.4 and 6 would lead to a stronger interleaflet en-
thalpic contribution for a lower position of unsaturation
due to the higher average S
CD
of those lipids.
We define the following term for the entropic behaviour
in the plane of a leaflet, where n
α,ik
is the count of lipids
of species α in the local region of site i in leaflet k shown
in fig.2(a), and φ
α,ik
is the fraction of lipids of species α
in the same local region.
H
entr
lat,i
=
V
S
k
B
T
4
n
α,ik
ln(φ
α,ik
) + n
α
0
,ik
ln(φ
α
0
,ik
)
(9)
This, summed over all sites of both leaflets, the above
equation can be expressed as:
H
entr
lat
=
V
S
k
B
T
4
X
k=0,1
X
i
n
α,ik
ln(φ
α,ik
)
+n
α
0
,ik
ln(φ
α
0
,ik
)
(10)
For the entropic interaction coupling the two leaflets,
we define
H
entr
inter,i
= V
0
S
k
B
T
1 +
P
p
α
l
α
n
0
α,ik
ln(φ
0
α,ik
)
+
1 +
P
p
α
0
l
α
0
n
0
α
0
,ik
ln(φ
0
α
0
,ik
)
(11)
.CC-BY 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted September 30, 2021. ; https://doi.org/10.1101/2021.09.29.462263doi: bioRxiv preprint
4
(a) (b)
FIG. 2: Cartoon representation of a local region of the
model, where the red sites are the locality over which
the parameters n
α
and φ
α
are calculated for a lattice
site i. (a) represents the sites for the entropic term of
H
lateral
, and (b) represents the sites for the entropic
term of H
interleaflet
which, when summed over all sites of both leaflets, be-
comes
H
entr
inter
= V
0
S
k
B
T
X
k=0,1
X
i
1 +
P
p
α
l
α
n
0
α,ik
ln(φ
0
α,ik
)
+
1 +
P
p
α
0
l
α
0
n
0
α
0
,ik
ln(φ
0
α
0
,ik
)
(12)
where p
α
/l
α
gives the fractional position of unsatura-
tion along the lipid tail, being 0 at the head group, and
1 at the end of the tail, which is then summed over all
unsaturations of the tail. This scaling of the entropic
terms ensures that the entropic contribution is stronger
for terms involving lipids with lower position of unsatu-
ration while the saturated lipids remain unaffected. n
0
α
and φ
0
α
are calculated in a locality shown in fig.2(b).
For a system of two lipids A and B, this gives us a
Hamiltonian with 9 parameters that are tunable during
simulations, shown together for convenience in table.I be-
low.
We use the Hamiltonian as defined above to calculate
the energy for our systems in the Monte-Carlo simula-
tions used in this study, which we describe in the follow-
ing section.
III. METHODS
We begin this section by describing how we obtained
the site variables and describing the systems we stud-
ied, followed by how they were used to generate data
using Monte-Carlo simulations. We also provide details
of the analysis tools used to obtain quantitative results
on phase separation and domain registration tendencies
of the systems from the generated data.
TABLE I: Description of all the tunable independent
parameters of the Hamiltonian for a 2-component
system with lipid species A and B.
Parameter Description
T Temperature
V
S
Lateral entropic term strength constant
V
0
S
Interleaflet entropic term strength constant
V
AA
Lateral enthalpic strength constant for A-A
interaction
V
BB
Lateral enthalpic strength constant for B-B
interaction
V
AB
Lateral enthalpic strength constant for A-B
interaction
V
0
AA
Interleaflet enthalpic strength constant for
A-A interaction
V
0
BB
Interleaflet enthalpic strength constant for B-
B interaction
V
0
AB
Interleaflet enthalpic strength constant for
A-B interaction
A. All-atom simulations of artificial lipids
We train the values of site variable (see eqn.8) for
each of the lipids used in the Monte-Carlo simulations
from all atom trajectories of DPPC-D23-Chol and
DPPC-D34-Chol systems. The study by S. Zhang and
X. Lin[4] used CGMD simulations where D23 and D34
are artificial lipids. To capture the lipid molecular prop-
erties faithfully, we reconstructed an atomic resolution
system and carried out AAMD simulations. In our
work, we have used a prescription that allows us to build
all-atom structures for artificial lipids that have a similar
chemical makeup to lipids with well-known forcefield
parameters. We provide a schematic of the workflow in
fig.3 and we describe the steps in the following text. We
prepared the all-atom descriptions of the artificial lipids,
namely D23 and D34 lipids using the Ligand Reader
and Modeler Module in CHARMM-GUI[33, 34]. Since
we aren’t using partial charges from Quantum-Chemical
calculations, we need to check the assigned partial
charges with a closely matching parameterized lipid
molecule which, in this case, were the partial charges
on the sp
2
and sp
3
hybridised lipid tail carbons in
the lipid DUPC. We then build the multi-component
lipid bilayer using the MemGen webserver[35] based
on the percentage composition of the required bilayer
and the number of lipids per leaflet in our system.
We have organized a tutorial with examples for the
workflow described above and made it publicly avail-
able at https://github.com/codesrivastavalab/
Membrane-Nanodomain-Registration/tree/main/
Artificial%20Lipid. Also, please note that common
modelling errors and issues during the generation
of the all-atom bilayer systems might require some
troubleshooting, which we discuss in sec.2 of the S.I.
For our systems, the average S
CD
was calculated for
.CC-BY 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted September 30, 2021. ; https://doi.org/10.1101/2021.09.29.462263doi: bioRxiv preprint
5
FIG. 3: Schematic of the workflow for generating
artificial lipid force field parameters for lipids with
different position of unsaturation and generation of the
bilayer topology.
DPPC, D23 and D34 lipids from the all-atom system of
144 DPPC, 80 D23/D34, and 56 Chol molecules in each
leaflet, which is close to the ratio used by Zhang and
Lin[4] in their publication. We used the CHARMM36
force field and first minimized the energy of the system
using steepest descent minimization. Then, we carried
out two rounds of 125 ps NVT equilibration at 303.15
K and four rounds of 250 ps NPT equilibration. We
ran production runs using GROMACS[36, 37] for a pe-
riod of 1µs with the PME[38, 39] method for electrostat-
ics, Nose-Hoover[40] temperature coupling and Parinello-
Rahman[41] pressure coupling. The LINCS[42] algorithm
was used to restrain hydrogen bonds. The tail order pa-
rameters of D23 and D34 lipids were calculated using the
gmx order program in GROMACS. The S
CD
values ob-
tained are averaged over all the D23/D34 lipids in the
D23/D34, DPPC, Chol ternary mixture and across 7.5
ns (S
CD
remain the same with increased averaging time).
S
CD
of n
th
carbon is calculated using the position of n−1
and n + 1 carbons.
Fig.4 shows the S
CD
values for D23 and D34 lipids.
The decrease in the S
CD
values occurs at the position of
double bonds in D23 and D34 lipids. The dip in value of
S
CD
for the C atoms near the end of the alkyl tails of D34
lipids indicates the lower position of the double bonds as
seen from the molecular structure of the lipids. On an
average, D34 lipids show more ordering of the alkyl tails
in comparison to D23 lipids. The final values of the site
variables calculated by averaging the S
CD
over all tail
carbon atoms for the DPPC-D23 systems were DPPC-
0.3074681, D23-0.1692048, and for the DPPC-D34 sys-
(a)
(b)
FIG. 4: a) shows the S
CD
for the tail carbon atoms of
the D23 and DPPC lipids as calculated by averaging
over all D23 lipids of the D23-DPPC system. b) Shows
the same for the D34-DPPC system.
tems, they were DPPC-0.3565843, D34-0.2400412. As
we can see, the average S
CD
is higher for the unsatu-
rated lipid with a lower position of unsaturation. Please
note that the purpose of the AAMD simulations is only
to extract out a more realistic site variable for our Hamil-
tonian.
B. Setup for Monte Carlo simulations
We studied 6 systems, each being a stack of 2 100 ×
100 square lattices. The upper and lower leaflet have
symmetric populations in each system, with no lattice
points left empty. So, each leaflet contains 10,000 lipids
.CC-BY 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted September 30, 2021. ; https://doi.org/10.1101/2021.09.29.462263doi: bioRxiv preprint
