Particlebased membrane model for mesoscopic simulation of cellular dynamics
TL;DR: In this article, a particlebased reactiondiffusion model for simulating lipid bilayer membranes is presented, based on interacting and reacting particles, each representing a coarse patch of a lipid monolayer and parameterized so as to reproduce the local membrane mechanics given by the Helfrich energy density over a range of relevant curvatures.
Abstract: We present a simple and computationally efficient coarsegrained and solventfree model for simulating lipid bilayer membranes. In order to be used in concert with particlebased reactiondiffusion simulations, the model is purely based on interacting and reacting particles, each representing a coarse patch of a lipid monolayer. Particle interactions include nearestneighbor bondstretching and anglebending, and are parameterized so as to reproduce the local membrane mechanics given by the Helfrich energy density over a range of relevant curvatures. Inplane fluidity is implemented with Monte Carlo bondflipping 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 particlebased 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 twodimensional shear flow under a gravity force is employed to measure the effective inplane viscosity of the membrane model, and show the possibility of modeling 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...
Summary (4 min read)
Jump to: [Introduction] – [II. THE MODEL] – [B. Differential geometry of the particlebased membrane model] – [C. Parameterspace optimization of interaction potentials] – [D. Additional interactions] – [E. Bondflipping moves] – [B. Time integration] – [C. Simulation code and visualization] – [A. Thermal undulations] – [B. Area compressibility] – [D. Nanoparticle wrapping] and [V. CONCLUSION]
Introduction
 The authors expect their model to be of high practical usability for ultra coarsegrained molecular dynamics or particlebased reactiondiffusion simulations of biological systems.
 Bilayer membranes have been the subject of computer simulations for more than three decades.33–35.
 Their approach consists of calibrating the interaction potentials to result in desired macroscopic mechanics and structural properties.
 The proposed model is essentially an elastic membrane model, comparable to triangulated models, with the difference that the desired elastic properties are reproduced through simple bonded interactions in contrast to complicated orientation or curvaturedependent potentials.
II. THE MODEL
 As shown in Fig. 1(a), two closepacked lattices of particles correspondingly represent the two leaflets of the membrane in this model.
 The elastic energy density contributed to the membrane is usually expressed in terms of the local curvature of the midsurface of the bilayer.
 The authors aim to avoid computing complex potential functions based on numerically obtained local curvature values.
 Thus, only bondstretching and anglebending 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.
B. Differential geometry of the particlebased membrane model
 In which the membrane is effectively composed of “particle dimers,” i.e., pairs of particles belonging to the top and bottom leaflets, a hypothetical midsurface is assumed to lie halfway between the particle dimers.
 Inspired by classical continuum shell theories, the authors assume that bending of the double layer deforms it such that a normal vector originating from a point p on the midsurface, pointing to a particle P on the upper or lower layer, remains perpendicular to the midsurface, independent of the state of deformation [see Fig. 1(b)].
 Without loss of generality, the authors focus on particles positioned on the top leaflet for the following derivations.
 For the purpose of calculating partial derivatives of the normal vector, Weingarten’s formula, n,µ = −bνµeν = −bµνgνσeσ has been used.
 Up to this point, the derived equations hold in all local coordinate systems at point p. A smart choice of the coordinate system can simplify the equations considerably.
C. Parameterspace optimization of interaction potentials
 Now that the authors have obtained equations describing the relative configuration of model particles in an arbitrarily curved membrane [Eqs. (4) and (5)], they can select interaction potentials which are functions of rPQ  and θpPQ and calculate effective energy densities corresponding to arbitrary curvature states.
 With such a definition for area elements, half of each lateral bond emanating from a particle P, plus all the outofplane angles having it as the vertex, is included in one area element around particle P.
 But without performing the simulation, the authors do not have a priori knowledge of the inplane angle, ψ, that this starshaped construct around each particle makes with the principal directions of the curvature of the midsurface.
 Thus, in general, the calculated effective energy density depends on this inplane angle.
 This way, the effective potential energy density is defined as feff = 〈 1 2 ∑ Ustretch + ∑ Ubending ∆A 〉 ψ , (10) where the summations are carried out for all interactions corresponding to one particle and ∆A denotes the area element.
D. Additional interactions
 The bondstretching and anglebending interactions described in Sec. II C only serve to reproduce the desired curvature elasticity of the membrane.
 In addition to these interactions, 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.
 Addition of this potential is necessary to hold the two leaflets together.
 Also, with the set of interactions described so far, volume exclusion is only present between neighboring membrane particles, and in principle, nonneighboring particles can interpenetrate.
E. Bondflipping moves
 In contrast, lipid bilayer membranes are twodimensional fluids in which lipid molecules can freely diffuse laterally, and this fluidity is essential for membrane remodeling.
 84 Following a scheme commonly used in triangulated membrane models,70,71 the inplane fluidity is introduced to the model via bondflipping Monte Carlo moves.
 This proposed move is accepted with the MetropolisHastings probability of exp [−β (Enew − Eold)] where Eold and Enew are the corresponding potential energies of the system in the old and new topologies, and β = 1/kT with k being the Boltzmann constant and T the temperature.
 In a simulation in the canonical ensemble, this lost energy will be compensated by the thermostat, which is the same as extracting work and adding equal amount of heat to the system.
 The frequency, φ, with which the bondflipping moves are proposed, acts as a control parameter for the model.
B. Time integration
 In order to simulate tensionless membranes in thermal equilibrium, an extended system dynamics approach is used to derive equations of motion and devise the proper numerical integration scheme.
 92–94 Thermostatting is achieved through NoséHoover chains, and isotropic cell fluctuations are used for barostatting to achieve zero inplane tension.
 On the other hand, in the absence of any solvent effects, and with the deterministic MTK integrator used here, the outofplane dynamics is solely determined by the particle masses and the stiffness of the forcefield developed based on the scheme introduced in Sec. II C.
 As the forcefield is the outcome of the parameterspace optimization aiming to reproduce the desired membrane elasticity, the only remaining parameter is the mass of model particles.
 Achieving physically relevant outofplane dynamics pertaining to membrane patches suspended in a solvent is only possible through either implementing a suitable stochastic integrator or including solvent effects explicitly or implicitly.
C. Simulation code and visualization
 Mainly due to the fact that the implementation of bondflipping Monte Carlo moves in available molecular dynamics software packages proved impractical, an inhouse C++ code has been developed to handle the simulations.
 Visualization is done via the Visual Molecular Dynamics (VMD) software package.95.
A. Thermal undulations
 A lipid bilayer patch in thermal equilibrium undergoes significant outofplane thermal undulations.
 These undulations can be studied from a statistical mechanics point of view to obtain energy distribution among different vibration modes.
 To ensure that the membrane patches have indeed been equilibrated, an estimate of the relaxation time of the system is required.
 The parameters for these fits are given in Table II.
 It is observed that increasing the lattice parameter in general has little effect on the ability of the model to reproduce continuum behavior.
B. Area compressibility
 As an example of additional physical properties of the model that can be taken into account when choosing potential parameters, area compressibility of the membrane is calculated for a model for which potential parameters are chosen from different families.
 The results are ordered as functions of the stiffness of the anglebending potential, Kb, and are depicted in Fig.
 The fluidity of the 2D liquid is described in terms of the surface viscosity, which arises from the assumed linear relation between the inplane shear stress and the corresponding velocity gradient.
 Thus, the general procedure described in this section has to be repeated if another integrator is used.
D. Nanoparticle wrapping
 As a final test of the usefulness of their membrane model to handle substantial deformations and model biologically relevant membrane remodeling processes, the authors simulate the interaction of spherical nanoparticles with the membrane, as a wellknown benchmark system.
 It is a useful test for the membrane model to show that (a) the model offers enough flexibility to simulate the budding behavior of bilayer membranes, and (b) if it correctly reproduces the interplay between bending and adhesion energies.
 For this type of nanoparticlemembrane interaction, and with a continuum membrane model, semianalytical studies105 have been carried out on the degree to which the surface of the nanoparticle is covered by the membrane, as a function of the dimensionless adhesion energy u = UpR2/κ as well as the potential range, ρ.
 It is observed that the model follows this prediction with very good accuracy.
 The good fit to the catenary curve is an indication that the particlebased model very well captures the zeroenergy regions and assumes corresponding minimal surface geometries.
V. CONCLUSION
 The authors have described a strongly coarsegrained model for simulating lipid bilayer membranes that is similar in nature with triangulated surface models, but is purely particle based, and as such is suitable for seamless integration into interactingparticle reactiondiffusion simulations.
 The model relies on bondstretching and anglebending interactions among nearestneighbor particles with parameters optimized to reproduce prescribed macroscopic curvature elasticity.
 These computer experiments have proven the model to be reliable in different equilibrium and nonequilibrium simulations and correctly predict the expected behavior of lipid bilayer membranes as twodimensional fluids obeying curvature elasticity.
 Armed with these capabilities, the authors ultimately aim to use this coarsegrained model in the context of iPRD simulations to study cellular signal transduction at large spatiotemporal scales.
Did you find this useful? Give us your feedback
Particlebased 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 (MesMP) model
The Journal of Chemical Physics 147, 044101 (2017); 10.1063/1.4993514
Enhanced configurational sampling with hybrid nonequilibrium molecular dynamics–Monte Carlo propagator
The Journal of Chemical Physics 148, 014101 (2018); 10.1063/1.5004154
Quantum mechanics/coarsegrained molecular mechanics (QM/CGMM)
The Journal of Chemical Physics 148, 014102 (2018); 10.1063/1.5006810
THE JOURNAL OF CHEMICAL PHYSICS 148, 044901 (2018)
Particlebased 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 BioSystems, 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 efﬁcient coarsegrained and solventfree model for simu
lating lipid bilayer membranes. In order to be used in concert with particlebased reactiondiffusion
simulations, the model is purely based on interacting and reacting particles, each representing a
coarse patch of a lipid monolayer. Particle interactions include nearestneighbor bondstretching and
anglebending and are parameterized so as to reproduce the local membrane mechanics given by the
Helfrich energy density over a range of relevant curvatures. Inplane ﬂuidity is implemented with
Monte Carlo bondﬂipping moves. The physical accuracy of the model is veriﬁed by ﬁve tests: (i)
Power spectrum analysis of equilibrium thermal undulations is used to verify that the particlebased
representation correctly captures the dynamics predicted by the continuum model of ﬂuid membranes.
(ii) It is veriﬁed 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 twodimensional shear ﬂow under a gravity force is employed
to measure the effective inplane viscosity of the membrane model and show the possibility of mod
eling membranes with speciﬁed 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 coarsegrained molecular
dynamics or particlebased reactiondiffusion 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@
fuberlin.de and frank.noe@fuberlin.de
b)
Electronic mail: thomas.weikl@mpikg.mpg.de
the simple case of equilibrating micronsized biomembranes,
a blind scale up in allatom molecular dynamics would be
out of reach of computational power for decades to come.
22
To ﬁll this computational gap and gain insights into cellular
processes, the development and application of coarsegrained
models are an important aspect of computer simulation. A
particularly promising framework to model cellular signaling
processes at membranes, involving space exclusions and spe
ciﬁc geometries found at membrane scaffolds, is particlebased
reactiondiffusion (PBRD) simulation,
21–24
especially the so
called interactingparticle reactiondiffusion (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 particlebased reactiondiffusion framework, espe
cially when it is required that the model be easily tunable,
robust, and yet computationally efﬁcient.
00219606/2018/148(4)/044901/11/$30.00 148, 0449011 Published by AIP Publishing.
0449012 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 speciﬁcally
developed forceﬁelds,
37
a vast variety of coarsegrained 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 coarsegraining achieved. This way, it becomes clear where
the proposed model ﬁts, and how it provides features suitable
for its integration in iPRD simulations.
The ﬁrst level of coarsegraining is achieved through
grouping a speciﬁc set of atoms in lipid molecules into
interaction centers and building effective forceﬁelds.
39
The
wellknown MARTINI forceﬁeld 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 coarsegrained mod
els for membranes composed of speciﬁc 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 curvatureinducing 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 coarsegraining pertains to oneparticlethick 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 sideeffect 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 particleparticle interactions.
59
Ayton
and Voth developed a systematic approach for parameteriz
ing the interaction potential in their EMDPD, 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, oneparticlethick models are also approached
as discretized continuum models. Triangulatedsurface 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 orientationdependent 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 ﬂuctuate, it is common practice with triangulatedsurface
models to utilize additional areapreserving constraints to con
trol the surface area of the membrane. Another approach is to
include the elasticity of an underlying continuous membrane
into a particlebased description through potentials that depend
on local surface ﬁtting.
76,77
Finally, the continuum descrip
tion with curvature elasticity can also be solved numerically
through available ﬁnite 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 coarsegrained mem
brane model which employs a twoparticlethick description
of the bilayer membrane, with each particle effectively rep
resenting a patch of lipids on each leaﬂet. This is a mini
mal structure that allows for ﬂexibility in modeling interac
tions of biomolecules with the membranes. The model relies
on simple bondstretching and anglebending potentials in
a dynamically updated bonded network and, thus, provides
enhanced computational efﬁciency through the exclusion of
nonbonded 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 curvaturedependent
potentials. Through a parameterspace 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 largescale simulations
of cellular dynamics and to speciﬁcally 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,
inplane viscosity, and budding under the inﬂuence of external
forces.
II. THE MODEL
As shown in Fig. 1(a), two closepacked lattices of par
ticles correspondingly represent the two leaﬂets 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 midsurface of the bilayer. We aim to avoid comput
ing complex potential functions based on numerically obtained
local curvature values. Thus, only bondstretching and angle
bending interactions amongst nearest neighbor particles are
considered. Considering an arbitrarily curved membrane, and
based on its local surface geometry, relative conﬁguration
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 ﬂuid bilayer
membrane is expressed as
79–81
0449013 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 leaﬂets in red and cyan, respectively. (b) Local surface geom
etry of the midsurface in an arbitrary state of deformation (blue surface) with
a collection of particle dimers whose positions are dictated by the midsurface
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 saddlesplay modulus. H and G represent
the mean and Gaussian curvatures, respectively, which are
deﬁned 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. Diﬀerential geometry of the particlebased
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 leaﬂets, a hypothetical midsurface 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 midsurface, pointing to a
particle P on the upper or lower layer, remains perpendicular
to the midsurface, 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 midsurface 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 leaﬂets, respectively. Without loss of
generality, we focus on particles positioned on the top leaﬂet
for the following derivations. For two neighboring particles
P and Q, corresponding midsurface 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
midsurface in the vicinity of point p. For point q, this descrip
tion can be approximated through a secondorder 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 ﬁrst 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 ﬁrst and second fundamental forms. The remaining
parameters in Eq. (4) are deﬁned 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 ﬁrst
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 deﬁnition:
u
σ
= u
?σ
−
1
2
Γ
σ
µν
u
?µ
u
?ν
, (8)
0449014 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
= δ
µ
ν
, ﬁrst order length differen
tials as well as the ﬁrst and second fundamental forms remain
unchanged.
C. Parameterspace optimization
of interaction potentials
Now that we have obtained equations describing the rel
ative conﬁguration 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 speciﬁc set of interaction poten
tials, as a function of midsurface 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 leaﬂets are connected via lateral bonds. Also, an angle
bending potential is assumed to exist for outofplane rotations
of such bonds. These two bonded interactions are handled,
respectively, with the following Morsetype bondstretching
and harmonic anglebending 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 leaﬂet) and P
0
is the particle residing
on the bottom leaﬂet, forming a dimer with particle P. The
equilibrium angle is chosen to be π/2, which corresponds to
anglebending with respect to a ﬂat 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 midsurface.
In order to calculate the effective energy density, an area
element on the midsurface corresponding to a set of interac
tions has to be deﬁned. We propose Voronoi tessellation be
used to do so in a systematic way. In the simple case of a
hexagonal closepacked lattice of particles, Voronoi tessella
tion simply yields hexagons centered at particles’ projections
on the midsurface. 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ﬂipping
Monte Carlo moves that will be discussed in Sec. II E. With
such a deﬁnition for area elements, half of each lateral bond
emanating from a particle P, plus all the outofplane 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 inplane angle, ψ, that this star
shaped construct around each particle makes with the principal
directions of the curvature of the midsurface. Thus, in gen
eral, the calculated effective energy density depends on this
inplane 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 deﬁned 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 leaﬂet, 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
deﬁned 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 midsurface 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 bondstretching and anglebending 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
tocenter 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 leaﬂets 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 ﬂip
between the leaﬂets. 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, nonneighboring particles can inter
penetrate. In simulations where this issue arises, a nonbonded
interaction, such as harmonic repulsion, can additionally be
included.
E. Bondﬂipping moves
The membrane model developed so far is based on a ﬁxed
topology of bonded interactions and, thus, pertains to a two
dimensional solid. In contrast, lipid bilayer membranes are
twodimensional ﬂuids in which lipid molecules can freely
diffuse laterally, and this ﬂuidity is essential for membrane
Citations
More filters
01 Mar 2003
TL;DR: In this paper, the physical properties of bilayer membranes were investigated by using a simple and efficient computer model, where the amphiphilic molecules were modeled as short rigid trimers with finite range pair interactions between them.
Abstract: We use a simple and efficient computer model to investigate the physical properties of bilayer membranes. The amphiphilic molecules are modeled as short rigid trimers with finite range pair interactions between them. The pair potentials have been designed to mimic the hydrophobic interactions, and to allow the simulation of the membranes without the embedding solvent as if the membrane is in vacuum. We find that upon decreasing the area density of the molecules the membrane undergoes a solid–fluid phase transition, where in the fluid phase the molecules can diffuse within the membrane plane. The surface tension and the bending modulus of the fluid membranes are extracted from the analysis of the spectrum of thermal undulations. At low area densities we observe the formation of pores in the membrane through which molecules can diffuse from one layer to the other. The appearance of the pores is explained using a simple model relating it to the area dependence of the free energy.
14 citations
••
TL;DR: This review collates together and discusses the various mechanicsbased mesoscopic models for proteinmediated membrane deformation studies and provides an elaborate description of a mesoscopic model where the membrane is modeled as a triangulated sheet and proteins are represented as either nematics or filaments.
7 citations
••
TL;DR: In this paper , the authors present an extensible framework, Mem3DG, for modeling 3D mechanochemical dynamics of membranes based on discrete differential geometry (DDG) on triangulated meshes.
Abstract: Biomembranes adopt varying morphologies that are vital to cellular functions. Many studies use computational modeling to understand how various mechanochemical factors contribute to membrane shape transformations. Compared with approximationbased methods (e.g., finite element method [FEM]), the class of discrete mesh models offers greater flexibility to simulate complex physics and shapes in three dimensions; its formulation produces an efficient algorithm while maintaining coordinatefree geometric descriptions. However, ambiguities in geometric definitions in the discrete context have led to a lack of consensus on which discrete mesh model is theoretically and numerically optimal; a bijective relationship between the terms contributing to both the energy and forces from the discrete and smooth geometric theories remains to be established. We address this and present an extensible framework, Mem3DG, for modeling 3D mechanochemical dynamics of membranes based on discrete differential geometry (DDG) on triangulated meshes. The formalism of DDG resolves the inconsistency and provides a unifying perspective on how to relate the smooth and discrete energy and forces. To demonstrate, Mem3DG is used to model a sequence of examples with increasing mechanochemical complexity: recovering classical shape transformations such as 1) biconcave disk, dumbbell, and unduloid; and 2) spherical bud on spherical, flatpatch membrane; investigating how the coupling of membrane mechanics with protein mobility jointly affects phase and shape transformation. As highresolution 3D imaging of membrane ultrastructure becomes more readily available, we envision Mem3DG to be applied as an endtoend tool to simulate realistic cell geometry under userspecified mechanochemical conditions.
4 citations
••
09 Nov 2018
TL;DR: This work presents recent innovations in the acceleration of molecular dynamics in GPUs to simulate nonHamiltonian systems and shows the performance of measurepreserving geometric integrator in the canonical ensemble, that is, at constant temperature.
Abstract: Molecular dynamics simulation is currently the theoretical technique eligible to simulate a wide range of systems from soft condensed matter to biological systems. However, of the excellent results that the technique has arrogated, this approach remains computationally expensive, but with the emergence of the new supercomputing technologies bases on graphics processing units graphical processing unitsbased systems GPUs, the perspective has changed. The GPUs allow performing large and complex simulations at a significantly reduced time. In this work, we present recent innovations in the acceleration of molecular dynamics in GPUs to simulate nonHamiltonian systems. In particular, we show the performance of measurepreserving geometric integrator in the canonical ensemble, that is, at constant temperature. We provide a validation and performance evaluation of the code by calculating the thermodynamic properties of a LennardJones fluid. Our results are in excellent agreement with reported data reported from literature, which were calculated with CPUs. The scope and limitations for performing simulations of highthroughput MD under rigorous statistical thermodynamics in the canonical ensemble are discussed and analyzed.
2 citations
References
More filters
••
TL;DR: VMD is a molecular graphics program designed for the display and analysis of molecular assemblies, in particular biopolymers such as proteins and nucleic acids, which can simultaneously display any number of structures using a wide variety of rendering styles and coloring methods.
Abstract: VMD is a molecular graphics program designed for the display and analysis of molecular assemblies, in particular biopolymers such as proteins and nucleic acids. VMD can simultaneously display any number of structures using a wide variety of rendering styles and coloring methods. Molecules are displayed as one or more "representations," in which each representation embodies a particular rendering method and coloring scheme for a selected subset of atoms. The atoms displayed in each representation are chosen using an extensive atom selection syntax, which includes Boolean operators and regular expressions. VMD provides a complete graphical user interface for program control, as well as a text interface using the Tcl embeddable parser to allow for complex scripts with variable substitution, control loops, and function calls. Full session logging is supported, which produces a VMD command script for later playback. Highresolution raster images of displayed molecules may be produced by generating input scripts for use by a number of photorealistic imagerendering applications. VMD has also been expressly designed with the ability to animate molecular dynamics (MD) simulation trajectories, imported either from files or from a direct connection to a running MD simulation. VMD is the visualization component of MDScope, a set of tools for interactive problem solving in structural biology, which also includes the parallel MD program NAMD, and the MDCOMM software used to connect the visualization and simulation programs. VMD is written in C++, using an objectoriented design; the program, including source code and extensive documentation, is freely available via anonymous ftp and through the World Wide Web.
46,130 citations
••
TL;DR: A theory of the elasticity of lipid bilayers is proposed and it is argued that in the case of vesicles (= closed bilayer films) the only elasticity controlling nonspherical shapes is that of curvature.
4,853 citations
••
TL;DR: An improved and extended version of the coarse grained lipid model is presented, coined the MARTINI force field, based on the reproduction of partitioning free energies between polar and apolar phases of a large number of chemical compounds to reproduce the free energies of these chemical building blocks.
Abstract: We present an improved and extended version of our coarse grained lipid model. The new version, coined the MARTINI force field, is parametrized in a systematic way, based on the reproduction of partitioning free energies between polar and apolar phases of a large number of chemical compounds. To reproduce the free energies of these chemical building blocks, the number of possible interaction levels of the coarsegrained sites has increased compared to those of the previous model. Application of the new model to lipid bilayers shows an improved behavior in terms of the stress profile across the bilayer and the tendency to form pores. An extension of the force field now also allows the simulation of planar (ring) compounds, including sterols. Application to a bilayer/cholesterol system at various concentrations shows the typical cholesterol condensation effect similar to that observed in all atom representations.
4,580 citations
••
TL;DR: In this paper, a modularly invariant equations of motion are derived that generate the isothermalisobaric ensemble as their phase space averages, and the resulting methods are tested on two problems, a particle in a onedimensional periodic potential and a spherical model of C60 in the solid/fluid phase.
Abstract: Modularly invariant equations of motion are derived that generate the isothermal–isobaric ensemble as their phase space averages. Isotropic volume fluctuations and fully flexible simulation cells as well as a hybrid scheme that naturally combines the two motions are considered. The resulting methods are tested on two problems, a particle in a one‐dimensional periodic potential and a spherical model of C60 in the solid/fluid phase.
4,282 citations
••
TL;DR: The presented lipid FF is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains and is anticipated to be of utility for simulations of pure lipid systems as well as heterogeneous systems including membrane proteins.
Abstract: A significant modification to the additive allatom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gellike structures of bilayers well above the gel transition temperature, selected torsional, LennardJones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from highlevel ab initio calculations on small molecules to semiempirical QM studies on a 1,2dipalmitoylsnphosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: ...
3,489 citations
Related Papers (5)
Frequently Asked Questions (2)
Q2. What are the contributions mentioned in the paper "Particlebased membrane model for mesoscopic simulation of cellular dynamics" ?
In this paper, a particlebased model for simulating lipid bilayer membranes is presented, which is similar in nature with triangulated surface models, but is purely particle based, and as such is suitable for seamless integration into interactingparticle reactiondiffusion simulations.