Particle-based membrane model for mesoscopic simulation of cellular dynamics
Mohsen Sadeghi, Thomas R. Weikl, and Frank Noé
Citation: The Journal of Chemical Physics 148, 044901 (2018); doi: 10.1063/1.5009107
View online: https://doi.org/10.1063/1.5009107
View Table of Contents: http://aip.scitation.org/toc/jcp/148/4
Published by the American Institute of Physics
Articles you may be interested in
Molecular dynamics based enhanced sampling of collective variables with very large time steps
The Journal of Chemical Physics 148, 024106 (2018); 10.1063/1.4999447
SSAGES: Software Suite for Advanced General Ensemble Simulations
The Journal of Chemical Physics 148, 044104 (2018); 10.1063/1.5008853
Perspective: Dissipative particle dynamics
The Journal of Chemical Physics 146, 150901 (2017); 10.1063/1.4979514
The mesoscopic membrane with proteins (MesM-P) model
The Journal of Chemical Physics 147, 044101 (2017); 10.1063/1.4993514
Enhanced configurational sampling with hybrid non-equilibrium molecular dynamics–Monte Carlo propagator
The Journal of Chemical Physics 148, 014101 (2018); 10.1063/1.5004154
Quantum mechanics/coarse-grained molecular mechanics (QM/CG-MM)
The Journal of Chemical Physics 148, 014102 (2018); 10.1063/1.5006810
THE JOURNAL OF CHEMICAL PHYSICS 148, 044901 (2018)
Particle-based membrane model for mesoscopic simulation
of cellular dynamics
Mohsen Sadeghi,
1,a)
Thomas R. Weikl,
2,b)
and Frank No
´
e
1,a)
1
Department of Mathematics and Computer Science, Freie Universit
¨
at Berlin, Arnimallee 6,
14195 Berlin, Germany
2
Department of Theory and Bio-Systems, Max Planck Institute of Colloids and Interfaces,
Science Park Golm, 14424 Potsdam, Germany
(Received 12 October 2017; accepted 5 January 2018; published online 23 January 2018)
We present a simple and computationally efficient coarse-grained and solvent-free model for simu-
lating lipid bilayer membranes. In order to be used in concert with particle-based reaction-diffusion
simulations, the model is purely based on interacting and reacting particles, each representing a
coarse patch of a lipid monolayer. Particle interactions include nearest-neighbor bond-stretching and
angle-bending and are parameterized so as to reproduce the local membrane mechanics given by the
Helfrich energy density over a range of relevant curvatures. In-plane fluidity is implemented with
Monte Carlo bond-flipping moves. The physical accuracy of the model is verified by five tests: (i)
Power spectrum analysis of equilibrium thermal undulations is used to verify that the particle-based
representation correctly captures the dynamics predicted by the continuum model of fluid membranes.
(ii) It is verified that the input bending stiffness, against which the potential parameters are optimized,
is accurately recovered. (iii) Isothermal area compressibility modulus of the membrane is calculated
and is shown to be tunable to reproduce available values for different lipid bilayers, independent of
the bending rigidity. (iv) Simulation of two-dimensional shear flow under a gravity force is employed
to measure the effective in-plane viscosity of the membrane model and show the possibility of mod-
eling membranes with specified viscosities. (v) Interaction of the bilayer membrane with a spherical
nanoparticle is modeled as a test case for large membrane deformations and budding involved in
cellular processes such as endocytosis. The results are shown to coincide well with the predicted
behavior of continuum models, and the membrane model successfully mimics the expected budding
behavior. We expect our model to be of high practical usability for ultra coarse-grained molecular
dynamics or particle-based reaction-diffusion simulations of biological systems. Published by AIP
Publishing. https://doi.org/10.1063/1.5009107
I. INTRODUCTION
Lipid bilayer membranes are integral parts of the machin-
ery of living cells. Apart from the obvious role of providing
a mechanical and chemical barrier for the cell, they form the
boundary of nearly all the organelles inside the eukaryotic cells
and also take part in cellular functions such as signal transduc-
tion.
1
Biologically relevant processes at membranes, such as
protein recruitment and insertion, assembly of protein scaf-
fold at membranes, and membrane remodeling, often involve
spatial scales from tens to hundreds of nanometers, and time
scales from milliseconds to minutes.
1
As an example, consider
endocytosis and exocytosis at plasma membranes.
2
While all-
atom molecular simulations are extremely successful for the
study of individual macromolecules and small complexes
3–9
and can reach thermodynamics and kinetics at very long time
scales with the aid of enhanced sampling methods and Markov
state modeling,
7,10–18
they have severe limitations in terms
of system sizes that can be sampled exhaustively.
19
Even for
a)
Authors to whom correspondence should be addressed: mohsen.sadeghi@
fu-berlin.de and frank.noe@fu-berlin.de
b)
Electronic mail: thomas.weikl@mpikg.mpg.de
the simple case of equilibrating micron-sized biomembranes,
a blind scale up in all-atom molecular dynamics would be
out of reach of computational power for decades to come.
22
To fill this computational gap and gain insights into cellular
processes, the development and application of coarse-grained
models are an important aspect of computer simulation. A
particularly promising framework to model cellular signaling
processes at membranes, involving space exclusions and spe-
cific geometries found at membrane scaffolds, is particle-based
reaction-diffusion (PBRD) simulation,
21–24
especially the so-
called interacting-particle reaction-diffusion (iPRD) models
that include interaction forces between particles.
25–29
The par-
ticles in such models typically represent entire proteins, protein
domains, or metabolites and thus represent a spatial resolution
of a few nanometers. Despite the success of such models in
simulating cellular signal transduction processes,
28,30–32
these
approaches are missing membrane mechanics in order to be
able to model signaling at biomembranes. In spite of the exten-
sive research on membrane models, there is arguably no readily
usable model at the same scale that is suited to be integrated
into such a particle-based reaction-diffusion framework, espe-
cially when it is required that the model be easily tunable,
robust, and yet computationally efficient.
0021-9606/2018/148(4)/044901/11/$30.00 148, 044901-1 Published by AIP Publishing.
044901-2 Sadeghi, Weikl, and No
´
e J. Chem. Phys. 148, 044901 (2018)
Bilayer membranes have been the subject of computer
simulations for more than three decades.
33–35
Apart from all-
atom simulations based on general purpose
36
or specifically
developed force-fields,
37
a vast variety of coarse-grained com-
putational models developed for bilayer membranes exist (see
Refs. 33 and 38 and references therein for an overview). While
we do not aim to provide a comprehensive review of all coarse-
grained membrane models, it is useful to look at important
modeling approaches and categorize them based on the level
of coarse-graining achieved. This way, it becomes clear where
the proposed model fits, and how it provides features suitable
for its integration in iPRD simulations.
The first level of coarse-graining is achieved through
grouping a specific set of atoms in lipid molecules into
interaction centers and building effective force-fields.
39
The
well-known MARTINI force-field falls into this category.
40–42
Although considerably reducing the number of particles, this
approach is still only suitable for simulating relatively small
scale processes.
43
More recently, Srivastava and Voth devised a
general approach for developing similar coarse-grained mod-
els for membranes composed of specific lipid molecules or
lipid mixtures. Their approach consists of calibrating the inter-
action potentials to result in desired macroscopic mechan-
ics and structural properties.
44
Simunovic et al. employed
this model to successfully simulate membrane remodeling
by curvature-inducing proteins.
45–47
The next step in coarse-
graining is to develop “bead” models, in which various lipid
molecules are represented not by interaction groups pertain-
ing to their atomistic representations but by a small chain
of generic particles.
48–57
A major challenge with these mod-
els is the choice of interaction potentials.
50
A higher level
of coarse-graining pertains to one-particle-thick models in
which the curvature elasticity is recovered through orientation-
dependent pairwise interactions.
58–61
Drouffe et al. pioneered
this approach and showed that through these orientation-
dependent interactions, stable membranes and vesicles can
form,
58
although with the side-effect of predicting consider-
ably low bending rigidities. To control the bending rigidity,
Kohyama proposed a model in which local curvature of the
membrane affects the particle-particle interactions.
59
Ayton
and Voth developed a systematic approach for parameteriz-
ing the interaction potential in their EM-DPD, and later, EM2
membrane models.
62–65
They performed detailed atomistic
simulations and employed energy equivalence in bending and
bulk expansion/contraction modes to obtain optimal parame-
ters for the mesoscopic model. They further applied these mod-
els in the study of membrane remodeling.
65,66
From a different
perspective, one-particle-thick models are also approached
as discretized continuum models. Triangulated-surface mod-
els developed by Gompper and Kroll
67–69
and Noguchi and
Gompper
71,72
follow such an approach and, instead of rely-
ing on pairwise orientation-dependent potentials, use angle-
bending potentials between neighboring triangles to directly
reproduce the curvature elasticity in a discretized model.
Bahrami et al. used a similar model to study the interaction
of nanoparticles with membranes
72,73
and formation of mem-
brane tubules.
74
Atilgan and Sun also incorporated the effect
of transmembrane proteins into a triangulated model.
75
As the
dimensions and areas of triangular elements in these models
can fluctuate, it is common practice with triangulated-surface
models to utilize additional area-preserving constraints to con-
trol the surface area of the membrane. Another approach is to
include the elasticity of an underlying continuous membrane
into a particle-based description through potentials that depend
on local surface fitting.
76,77
Finally, the continuum descrip-
tion with curvature elasticity can also be solved numerically
through available finite element methods developed for thin
shell mechanics.
78
In effect, these approaches substitute par-
ticles with computational nodes of a discretized continuum
model.
In this paper, we introduce a novel coarse-grained mem-
brane model which employs a two-particle-thick description
of the bilayer membrane, with each particle effectively rep-
resenting a patch of lipids on each leaflet. This is a mini-
mal structure that allows for flexibility in modeling interac-
tions of biomolecules with the membranes. The model relies
on simple bond-stretching and angle-bending potentials in
a dynamically updated bonded network and, thus, provides
enhanced computational efficiency through the exclusion of
non-bonded pairwise interactions. The proposed model is
essentially an elastic membrane model, comparable to trian-
gulated models, with the difference that the desired elastic
properties are reproduced through simple bonded interactions
in contrast to complicated orientation- or curvature-dependent
potentials. Through a parameter-space optimization scheme,
these interactions are easily tuned to reproduce membranes
with desired elastic properties. The ultimate aim of develop-
ing such a model is to include it in large-scale simulations
of cellular dynamics and to specifically use it for studying
cellular signal transduction using iPRD models. The com-
puter experiments laid out in the following are designed to
show that, despite its relative simplicity, inexpensive simu-
lations done with the model very well reproduce expected
behavior in terms of thermal undulations, area compressibility,
in-plane viscosity, and budding under the influence of external
forces.
II. THE MODEL
As shown in Fig. 1(a), two close-packed lattices of par-
ticles correspondingly represent the two leaflets of the mem-
brane in this model. The elastic energy density contributed to
the membrane is usually expressed in terms of the local curva-
ture of the mid-surface of the bilayer. We aim to avoid comput-
ing complex potential functions based on numerically obtained
local curvature values. Thus, only bond-stretching and angle-
bending interactions amongst nearest neighbor particles are
considered. Considering an arbitrarily curved membrane, and
based on its local surface geometry, relative configuration
of particles and the resulting bond lengths and angles are
obtained. An effective energy density pertaining to bonded
interactions is thus calculated and compared with the cur-
vature elasticity modeled via the Helfrich energy density to
parameterize the interaction potentials.
A. Curvature elasticity of the bilayer membrane
The Helfrich energy density of a curved fluid bilayer
membrane is expressed as
79–81
044901-3 Sadeghi, Weikl, and No
´
e J. Chem. Phys. 148, 044901 (2018)
FIG. 1. (a) Snapshot of the proposed membrane model with particles forming
top and bottom leaflets in red and cyan, respectively. (b) Local surface geom-
etry of the mid-surface in an arbitrary state of deformation (blue surface) with
a collection of particle dimers whose positions are dictated by the mid-surface
geometry. Distances and angles between these particles are used in order to
probe the local curvature and relate between the particle model and continuum
description of membrane mechanics.
f
H
= 2κ(H − H
0
)
2
+ ¯κG, (1)
in which the constant κ is the bending rigidity or splay
modulus of the membrane and ¯κ is its Gaussian curva-
ture rigidity or saddle-splay modulus. H and G represent
the mean and Gaussian curvatures, respectively, which are
defined based on the principal curvatures, c
1
and c
2
, as
H =
(
c
1
+ c
2
)
/2 and G = c
1
c
2
. H
0
is the spontaneous mean
curvature of the membrane, corresponding to a local curva-
ture that is induced in the membrane not by external forces
but by internal effects such as the geometry of phospholipid
molecules.
82
B. Differential geometry of the particle-based
membrane model
In this model, in which the membrane is effectively com-
posed of “particle dimers,” i.e., pairs of particles belonging
to the top and bottom leaflets, a hypothetical mid-surface is
assumed to lie halfway between the particle dimers. Inspired
by classical continuum shell theories, we assume that bend-
ing of the double layer deforms it such that a normal vector
originating from a point p on the mid-surface, pointing to a
particle P on the upper or lower layer, remains perpendicular
to the mid-surface, independent of the state of deformation
[see Fig. 1(b)]. Thus, the position of the particle P is always
given as r
P
= r
p
±
t
2
n, where n is the normal vector of
the mid-surface at point p, t is the thickness of the mem-
brane, and the plus and minus signs correspond to particles
on the top or bottom leaflets, respectively. Without loss of
generality, we focus on particles positioned on the top leaflet
for the following derivations. For two neighboring particles
P and Q, corresponding mid-surface projections are consid-
ered to be p and q, given by local coordinates, r
p
=
(
0, 0
)
and r
q
=
u
1
, u
2
, respectively [Fig. 1(b)]. Thus, the set of
coordinates, u
1
and u
2
, provide a local parameterization of the
mid-surface in the vicinity of point p. For point q, this descrip-
tion can be approximated through a second-order Taylor
expansion
r
q
≈ r
p
+ u
µ
e
µ
+
1
2
u
µ
u
ν
e
µ,ν
= r
p
+ u
µ
e
µ
+
1
2
u
µ
u
ν
Γ
σ
µν
e
σ
+ b
µν
n
, (2)
where e
µ
= ∂
µ
r are the base vectors for the tangent space
at point p, Γ
σ
µν
= e
µ,ν
· e
σ
are the Christoffel symbols of
the second kind, and b
µν
= e
µ,ν
·n are the components of the
second fundamental form tensor.
83
It is to be noted that sum-
mation convention between a pair of upper and lower indices
is used here. Similarly, another Taylor expansion can be used
to approximate the normal vector at point q, making it possible
to express the position of particle Q with respect to particle P
as
r
PQ
= r
Q
− r
P
≈
"
u
σ
+
1
2
Γ
σ
µν
u
µ
u
ν
−
h
2
b
µν
g
νσ
u
µ
#
e
σ
+
1
2
b
µν
u
µ
u
ν
n, (3)
in which g
µν
= e
µ
·e
ν
is the metric tensor and we have g
µσ
g
σν
= g
µσ
g
σν
= δ
µ
ν
with δ
µ
ν
being Kronecker’s delta. For the
purpose of calculating partial derivatives of the normal vec-
tor, Weingarten’s formula, n
,µ
= −b
ν
µ
e
ν
= −b
µν
g
νσ
e
σ
has
been used.
83
It is noteworthy to mention that first order partial
derivatives of the normal vector contain second order deriva-
tives of the position vector, r, through the inclusion of the
b
µν
tensor components, effectively making the two approxi-
mations of the same order. The length of the vector r
PQ
as
well as the angle it makes with the normal vector at point p are
obtained by forming the following inner products:
|r
PQ
|
2
= r
PQ
· r
PQ
≈ I(u)
1 −
t
2
4
G
!
− II(u) t
1 −
t
2
H
+
1
4
II
2
(u) + C
1
−
t
2
C
2
+
1
4
C
3
(4)
and
|r
PQ
|
|
n
|
cos θ
p PQ
= −r
PQ
· n ≈ −
1
2
II(u), (5)
where
I(u) = g
µν
u
µ
u
ν
II(u) = b
µν
u
µ
u
ν
(6)
are the first and second fundamental forms. The remaining
parameters in Eq. (4) are defined as follows:
C
1
= Γ
µνσ
u
µ
u
ν
u
σ
,
C
2
= Γ
σ
µν
b
σγ
u
µ
u
ν
u
γ
,
C
3
= Γ
σ
µν
Γ
µ
0
ν
0
σ
u
µ
u
ν
u
µ
0
u
ν
0
, (7)
where Γ
µνσ
= Γ
γ
µν
g
γσ
are the Christoffel symbols of the first
kind.
83
Up to this point, the derived equations hold in all local
coordinate systems at point p. A smart choice of the coordi-
nate system can simplify the equations considerably. A perfect
candidate is the locally tangent coordinate system, with the
following implicit definition:
u
σ
= u
?σ
−
1
2
Γ
σ
µν
u
?µ
u
?ν
, (8)
044901-4 Sadeghi, Weikl, and No
´
e J. Chem. Phys. 148, 044901 (2018)
in which
u
?1
, u
?2
are the new coordinates with the same ori-
gin at point p, and the Christoffel symbols are calculated in the
old coordinate system at point p. It can be shown that in this
coordinate system, Christoffel symbols vanish identically, and
yet, because
∂u
µ
/∂u
?ν
p
= δ
µ
ν
, first order length differen-
tials as well as the first and second fundamental forms remain
unchanged.
C. Parameter-space optimization
of interaction potentials
Now that we have obtained equations describing the rel-
ative configuration of model particles in an arbitrarily curved
membrane [Eqs. (4) and (5)], we can select interaction poten-
tials which are functions of |r
PQ
| and θ
pPQ
and calculate
effective energy densities corresponding to arbitrary curva-
ture states. In effect, we seek to obtain numerical values of the
energy density arising from a specific set of interaction poten-
tials, as a function of mid-surface curvature, prior to running
an actual simulation with these potentials. As a simple choice,
we assume that nearest neighbor particles on both top and
bottom leaflets are connected via lateral bonds. Also, an angle-
bending potential is assumed to exist for out-of-plane rotations
of such bonds. These two bonded interactions are handled,
respectively, with the following Morse-type bond-stretching
and harmonic angle-bending potentials:
U
stretch
|r
PQ
|
= D
e
f
1 − e
−α
(
|
r
PQ
|
−a
)
g
2
U
bend
θ
P
0
P Q
= K
b
θ
P
0
PQ
−
π
2
2
, (9)
where a denotes the lattice parameter (or equilibrium separa-
tion of particles on each leaflet) and P
0
is the particle residing
on the bottom leaflet, forming a dimer with particle P. The
equilibrium angle is chosen to be π/2, which corresponds to
angle-bending with respect to a flat membrane.
It is to be noted that this choice of bonded interactions,
and the potentials to handle them, is by no means unique. The
general procedure laid out here can be applied to many other
choices, with the condition that geometric information can be
extracted uniquely from the curvature of the mid-surface.
In order to calculate the effective energy density, an area
element on the mid-surface corresponding to a set of interac-
tions has to be defined. We propose Voronoi tessellation be
used to do so in a systematic way. In the simple case of a
hexagonal close-packed lattice of particles, Voronoi tessella-
tion simply yields hexagons centered at particles’ projections
on the mid-surface. Although in general, the shape and area
of elements corresponding to particle projections are a func-
tion of their coordination number, especially considering the
fact that the coordination number changes due to bond-flipping
Monte Carlo moves that will be discussed in Sec. II E. With
such a definition for area elements, half of each lateral bond
emanating from a particle P, plus all the out-of-plane angles
having it as the vertex, is included in one area element around
particle P. But without performing the simulation, we do not
have a priori knowledge of the in-plane angle, ψ, that this star-
shaped construct around each particle makes with the principal
directions of the curvature of the mid-surface. Thus, in gen-
eral, the calculated effective energy density depends on this
in-plane angle. To compensate for this ambiguity, and avoid
directional bias, the effective energy density is numerically
averaged out over possible values of ψ. This way, the effective
potential energy density is defined as
f
eff
=
*
1
2
P
U
stretch
+
P
U
bending
∆A
+
ψ
, (10)
where the summations are carried out for all interactions cor-
responding to one particle and ∆A denotes the area element.
The same procedure applies to the pair particle, P
0
, which lies
on the bottom leaflet, and the corresponding energy density is
simply added to f
eff
.
The chosen interaction potentials given in Eq. (9) contain
a set of parameters, D
e
, α, and K
b
. In order to obtain optimal
values for these parameters, a dimensionless error measure is
defined as
e =
∫
dc
1
∫
dc
2
(
f
eff
− f
H
)
2
∫
dc
1
∫
dc
2
f
2
H
, (11)
in which the integration is carried out in the mid-surface cur-
vature space spanned by its principal curvatures, c
1
and c
2
.
The integration range is arbitrary and corresponds to the range
of curvatures that have practical relevance. Minimizing this
error measure with respect to potential parameters yields their
optimal values.
D. Additional interactions
The bond-stretching and angle-bending interactions
described in Sec. II C only serve to reproduce the desired
curvature elasticity of the membrane. In addition to these inter-
actions, other potentials can be added to the model for different
geometrical or mechanical considerations, as long as they do
not perturb the effective energy density described by Eq. (10).
Most importantly, a harmonic potential is added between
particles in each dimer. This potential is described by the
expression U
thickness
(
t
)
= K
t
(
t − t
0
)
2
where t is the center-
to-center distance of the particles, and t
0
is the prescribed
equilibrium thickness of the membrane. Addition of this poten-
tial is necessary to hold the two leaflets together. The potential
strength, K
t
, should be chosen high enough to preserve the
thickness variations within a reasonable range and to pre-
vent the thermal motion of the particles to cause them to flip
between the leaflets. Note that this potential may also be cali-
brated with respect to the actual stiffness of bilayer membranes
across their thickness.
Also, with the set of interactions described so far, vol-
ume exclusion is only present between neighboring membrane
particles, and in principle, non-neighboring particles can inter-
penetrate. In simulations where this issue arises, a non-bonded
interaction, such as harmonic repulsion, can additionally be
included.
E. Bond-flipping moves
The membrane model developed so far is based on a fixed
topology of bonded interactions and, thus, pertains to a two-
dimensional solid. In contrast, lipid bilayer membranes are
two-dimensional fluids in which lipid molecules can freely
diffuse laterally, and this fluidity is essential for membrane