A stabilized finite element formulation for liquid shells and its
application to lipid bilayers
Roger A. Sauer∗1 , Thang X. Duong∗ , Kranthi K. Mandadapu†§ and David J. Steigmann‡
∗

†

Aachen Institute for Advanced Study in Computational Engineering Science (AICES),
RWTH Aachen University, Templergraben 55, 52056 Aachen, Germany

Department of Chemical and Biomolecular Engineering, University of California at Berkeley,
110 Gilman Hall, Berkeley, CA 94720-1460, USA
§

Chemical Sciences Division, Lawrence Berkeley National Laboratory, CA 94720, USA

arXiv:1601.03907v2 [] 18 Jan 2016

‡

Department of Mechanical Engineering, University of California at Berkeley,
6141 Etcheverry Hall, Berkeley, CA 94720-1740, USA

Abstract
This paper presents a new finite element (FE) formulation for liquid shells that is based on an
explicit, 3D surface discretization using C 1 -continuous finite elements constructed from NURBS
interpolation. Both displacement-based and mixed FE formulations are proposed. The latter is
needed for area-incompressible material behavior, where penalty-type regularizations can lead
to misleading results. In order to obtain quasi-static solutions, several numerical stabilization
schemes are proposed based on either stiffness, viscosity or projection. Several numerical examples are considered in order to illustrate the accuracy and the capabilities of the proposed
formulation, and to compare the different stabilization schemes. The presented formulation is
capable of simulating non-trivial surface shapes associated with tube formation and proteininduced budding of lipid bilayers. In the latter case, the presented formulation yields nonaxisymmetric solutions, which have not been observed in previous simulations. It is shown that
those non-axisymmetric shapes are preferred over axisymmetric ones.
Keywords: cell budding, cell tethering, Helfrich energy, isogeoemetric analysis, non-linear
finite elements, non-linear shell theory

Contents
1 Introduction

4

2 Summary of thin liquid shell theory

6

2.1

Thin shell kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2

Quasi-static equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3

Constitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3.1

Area-compressible lipid bilayer . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3.2

Area-incompressible lipid bilayer . . . . . . . . . . . . . . . . . . . . . . .

9

2.3.3

Model properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.4
1

Weak form

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

corresponding author, email: sauer@aices.rwth-aachen.de

1

3 Liquid shell stabilization
3.1

12

Adding stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.1.1

In-plane shear and bulk stabilization . . . . . . . . . . . . . . . . . . . . . 12

3.1.2

Sole in-plane shear stabilization . . . . . . . . . . . . . . . . . . . . . . . . 12

3.1.3

Relation to physical viscosity . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.2

Normal projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.3

Summary of the stabilization schemes . . . . . . . . . . . . . . . . . . . . . . . . 15

4 FE formulation

15

4.1

FE approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4.2

Discretized weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

4.3

Area constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.4

Solution procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.5

Rotational constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

4.6

Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5 Numerical examples

18

5.1

Pure bending and stretching of a flat strip . . . . . . . . . . . . . . . . . . . . . . 18

5.2

Inflation of a sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

5.3

Drawing of a tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

5.4

Cell budding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.4.1

Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5.4.2

Computational setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5.4.3

Bud shapes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

5.4.4

Surface energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5.4.5

Surface tension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5.4.6

Effective shear stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.4.7

Influence of the area-compressibility . . . . . . . . . . . . . . . . . . . . . 34

6 Conclusion

36

e
A Linearization of finto

37

2

B Linearization of ge

38

List of important symbols
1
aα
Aα
aα
Aα
aα,β
aα;β
aαβ
Aαβ
a
A
bαβ
Bαβ
B
cαβγδ
C
γ
Γγαβ
da
dA
δ...
e

fe
g
G
Ge
ge
H
H0
I
I1 , I2
i
I
J
k
k∗
K
Keff
ke
κ
κ1 , κ2
LI
λ1 , λ2
me

identity tensor in R3
co-variant tangent vectors of surface S at point x; α = 1, 2
co-variant tangent vectors of surface S0 at point X; α = 1, 2
contra-variant tangent vectors of surface S at point x; α = 1, 2
contra-variant tangent vectors of surface S0 at point X; α = 1, 2
parametric derivative of aα w.r.t. ξ β
co-variant derivative of aα w.r.t. ξ β
co-variant metric tensor components of surface S at point x
co-variant metric tensor components of surface S0 at point X
class of stabilization methods based on artificial shear viscosity
class of stabilization methods based on artificial shear stiffness
co-variant curvature tensor components of surface S at point x
co-variant curvature tensor components of surface S0 at point X
left surface Cauchy-Green tensor
contra-variant components of the material tangent
right surface Cauchy-Green tensor
surface tension of S
Christoffel symbols of the second kind
differential surface element on S
differential surface element on S0
variation of ...
index numbering the finite elements; e = 1, ..., nel
penalty parameter
finite element force vector of element Ωe
expression for the area-incompressibility constraint
expression for the weak form
contribution to G from finite element Ωe
finite element ‘force vector’ of element Ωe due to constraint g
mean curvature of S at x
spontaneous curvature prescribed at x
index numbering the finite element nodes
first and second invariants of the surface Cauchy Green tensors
surface identity tensor on S
surface identity tensor on S0
surface area change
bending modulus
Gaussian modulus
initial in-plane membrane bulk modulus
effective in-plane membrane bulk modulus
finite element tangent matrix associated with f e and ge
Gaussian curvature of surface S at x
principal curvatures of surface S at x
pressure shape function of finite element node I
principal surface stretches of S at x
number of pressure nodes of finite element Ωe
3

mν , mτ
m̄ν , m̄τ
M αβ
µ
µeff
nno
nel
nmo
ne
N αβ
NI
n
N
N
ν
ν
ξα
P
q
q
qe
Sα
S
S0
σ
σ αβ
t
t̄
T
Tα
V, Q
ϕ
ϕ̄
w
W
x
X
xI
XI
x
xe
Xe
Ωe
Ωe0

1

bending moment components acting at x ∈ ∂S
prescribed bending moment components
contra-variant bending moment components
initial in-plane membrane shear stiffness
effective in-plane membrane shear stiffness
total number of finite element nodes used to discretize S
total number of finite elements used to discretize S
total number of finite element nodes used to discretize pressure q
number of displacement nodes of finite element Ωe
total, contra-variant, in-plane membrane stress components
displacement shape function of finite element node I
surface normal of S at x
surface normal of S0 at X
array of the shape functions for element Ωe
in-plane membrane shear viscosity
normal vector on ∂S
convective surface coordinates; α = 1, 2
class of stabilization methods based on normal projection; projection matrix
Lagrange multiplier associated with area-incompressibility
array of all Lagrange multipliers qI in the system; I = 1, ..., nmo
array of all Lagrange multipliers qI for finite element Ωe ; I = 1, ..., me
contra-variant, out-of-plane shear stress components
current configuration of the surface
initial configuration of the surface
Cauchy stress tensor of the shell
stretch related, contra-variant, in-plane membrane stress components
effective traction acting on the boundary ∂S normal to ν
prescribed boundary tractions on Neumann boundary ∂t S
traction acting on the boundary ∂S normal to ν
traction acting on the boundary ∂S normal to aα
admissible function spaces
deformation map of surface S
prescribed boundary deformations on boundary ∂x S
hyperelastic stored surface energy density per current surface area
hyperelastic stored surface energy density per reference surface area
current position of a surface point on S
initial position of x on the reference surface S0
position vector of finite element node I lying on S
initial position of finite element node I on S0
array of all nodal positions xI of the discretized surface; I = 1, ..., nno
array of all nodal positions xI for finite element Ωe ; I = 1, ..., ne
array of all nodal positions X I for finite element Ωe0 ; I = 1, ..., ne
current configuration of finite element e
reference configuration of finite element e

Introduction

Biological membranes form the boundaries of cells and cell-internal organelles such as the endoplasmic reticulum, the golgi complex, mitochondria and endosomes. Mechanically they are
4

liquid shells that exhibit fluid-like behavior in-plane and solid-like behavior out-of-plane. They
mainly consist of self-assembled lipid bilayers and proteins. At the macroscopic level, these
membranes exist in different shapes such as invaginations, buds and cylindrical tubes (Zimmerberg and Kozlov, 2006; Shibata et al., 2009; McMahon and Gallop, 2005; Shemesh et al.,
2014). These shapes arise as a result of the lateral loading due to cytoskeletal filaments and
protein-driven spontaneous curvature. Cell membranes undergo many morphological and topological shape transitions to enable important biological processes such as endocytosis (Buser
and Drubin, 2013; Kukulski et al., 2012; Peter et al., 2004), cell motility (Keren, 2011) and
vesicle formation (Gruenberg and Stenmark, 2004; Budin et al., 2009). The shape transitions
occur as a result of lateral loading on the membranes from cytoskeletal filaments, such as actin,
from osmotic pressure gradients across the membrane, or from membrane-protein interactions.
For example, in endocytosis – a primary mode of transport of cargo between the exterior of
the cell and its interior – proteins bind to the flat membrane and induce invaginations or bud
shapes, followed by a tube formation by actin-mediated pulling forces (Kukulski et al., 2012;
Walani et al., 2015).
Most of the computational studies regarding biological membranes are restricted to shapes
resulting from axisymmetric conditions. However, many of the processes in cells are nonaxisymmetric in nature. To the best of our knowledge, there exist only a few studies that allow
for general, non-axisymmetric shapes. Therefore, it is important to advance computational
methods that can yield solutions for general conditions.
In the past, several computational models have been proposed for cell membranes. Depending
on how the membrane is discretized, two categories can be distinguished: Models based on
an explicit surface discretization, and models based on an implicit surface discretization. In
the second category, the surface is captured by a phase field (Du and Wang, 2007) or level
set function (Salac and Miksis, 2011) that is defined on the surrounding volume mesh. In the
first category, the surface is captured directly by a surface mesh. The approach is particularly
suitable if only surface effects are studied, such that no surrounding volume mesh is needed.
This is the approach taken here. An example is to use Galerkin surface finite elements: The
first corresponding 3D FE model for lipid bilayer membranes seems to be the formulation
of Feng and Klug (2006) and Ma and Klug (2008). Their FE formulation is based on socalled subdivision surfaces (Cirak and Ortiz, 2001), which provide C 1 -continuous FE surface
discretizations. Such discretizations are advantageous, since they do not require additional
degrees of freedom as C 0 -continuous FE formulations do. Still, C 0 -continuous FEs have been
considered to model red blood cell (RBC) membranes and their supporting protein skeleton
(Dao et al., 2003; Peng et al., 2010), phase changes of lipid bilayers (Elliott and Stinner, 2010),
and viscous cell membranes (Tasso and Buscaglia, 2013). Subdivision finite elements have been
used to study confined cells (Kahraman et al., 2012). Lipid bilayers can also be modeled with socalled ‘solid shell’ (i.e. classical volume) elements instead of surface shell elements (Kloeppel and
Wall, 2011). Using solid elements, C 0 -continuity is sufficient, but the formulation is generally
less efficient. For two-dimensional and axisymmetric problems also C 1 -continuous B-Spline and
Hermite finite elements have been used to study membrane surface flows (Arroyo and DeSimone,
2009; Rahimi and Arroyo, 2012), cell invaginations (Rim et al., 2014), and cell tethering and
adhesion (Rangarajan and Gao, 2015). The latter work also discusses the generalization to
three-dimensional B-spline FE. For some problems it is also possible to use specific, Mongepatch FE discretizations (Rangamani et al., 2013, 2014).
There are also several works that do not use finite element approaches. Examples are numerical
ODE integration (Agrawal and Steigmann, 2009), Monte Carlo methods (Ramakrishnan et al.,
2010), molecular dynamics (Li and Lykotrafitis, 2012), finite difference methods (Lau et al.,
2012; Gu et al., 2014) and mesh-free methods (Rosolen et al., 2013). There are also non5

Galerkin FE approaches that use triangulated surfaces, e.g. see Jarić et al. (1995); Jie et al.
(1998).
For quasi-static simulations of liquid membranes and shells, the formulation needs to be stabilized. Therefore, various stabilization methods have been proposed considering artificial viscosity (Ma and Klug, 2008; Sauer, 2014), artificial stiffness (Kahraman et al., 2012) and normal
offsets – either as a projection of the solution (with intermediate mesh update steps) (Sauer,
2014), or as a restriction of the surface variation (Rangarajan and Gao, 2015). The instability
problem is absent, if shear stiffness is present, e.g. due to an underlying cytoskeleton, like in
RBCs (Dao et al., 2003; Peng et al., 2010; Kloeppel and Wall, 2011). However, this setting does
not apply to (purely) liquid membranes and shells.
Here, a novel FE formulation is presented for liquid shells that is based on an explicit, 3D
surface discretization using NURBS (Cottrell et al., 2009). The following aspects are new:
• A theoretical study of the effective bulk and shear stiffness of the Helfrich bending model,
• the use of 3D, C 1 -continuous NURBS-based surface discretizations, considering purely
displacement-based and
• mixed, LBB-conforming2 finite elements,
• a comparison of both these FE formulations, illustrating the limitations of the former,
• the formulation, application and comparison of various stabilization schemes for quasistatic simulations
• that allow to accurately compute the in-plane stress state of the liquid, like the surface
tension,
• the verification of the formulation by various analytical solutions,
• the possibility to capture complex, non-axisymmetric solutions, and
• new insight into the budding of cells and vesicles.
To the best of our knowledge, non-axisymmetric shapes have not been simulated to the same
detail as is presented here.
The remainder of this paper is organized as follows: Sec. 2 presents an overview of the underlying
theory of thin liquid shells. The governing weak form is presented and the model properties are
discussed. Stabilization methods are then addressed in Sec. 3. The finite element formulation
follows in Sec. 4. The formulation is verified and illustrated by several 3D numerical examples
in Sec. 5. In particular, tethering and budding are studied. The paper then concludes with
Sec. 6.

2

Summary of thin liquid shell theory

This section gives a brief overview of the governing equations of thin shell theory considering
liquid material behavior. A much more detailed discussion can be found in Sauer and Duong
(2015) and in the references therein.
2

satisfying the Ladyzhenskaia-Babuška-Brezzi condition, see e.g. Babuška (1973); Bathe (1996)

6

2.1

Thin shell kinematics

The surface of the shell, denoted by S, is described by the parameterization
x = x(ξ α ) ,

α = 1, 2,

(1)

where ξ 1 and ξ 2 are coordinates that can be associated with a flat 2D domain that is then
mapped to S by function (1). The surface is considered deformable. The deformation of
S is assessed in relation to a reference configuration S0 , described by the parameterization
X = X(ξ α ). S is considered here to coincide with S0 initially (e.g. at time t = 0), such that
x = X for t = 0. From these mappings, the tangent vectors Aα = X ,α and aα = x,α can be
determined. Here ...,α := ∂.../∂ξ α denotes the parametric derivative. From the tangent vectors,
the surface normals
A1 × A2
N=
,
(2)
kA1 × A2 k
a1 × a2
n=
,
(3)
ka1 × a2 k
and the metric tensor components
Aαβ = Aα · Aβ ,

(4)

aαβ = aα · aβ

(5)

can be defined. From the inverses [Aαβ ] := [Aαβ ]−1 and [aαβ ] := [aαβ ]−1 , the dual tangent
vectors
Aα = Aαβ Aβ ,
(6)
aα = aαβ aβ

(7)
α

(with summation implied on repeated indices) can be introduced, such that A · Aβ = δβα and
aα · aβ = δβα , where δβα is the Kronecker symbol. Further, we can define the surface identity
tensors
I = Aα ⊗ Aα ,
(8)
i = aα ⊗ aα ,

(9)

such that the usual, 3D identity becomes
1=I +N ⊗N =i+n⊗n .

(10)

From the second derivatives of Aα and aα one can define the curvature tensor components
Bαβ = Aα,β · N ,

(11)

bαβ = aα,β · n ,

(12)

allowing us to compute the mean, Gaussian and principal curvatures
H = 21 bαβ aαβ ,

κ=

det[bαβ ]
,
det[aαβ ]

κ1/2 = H ±

p
H2 − κ

(13)

of surface S, and likewise for S0 . Here, bαβ = aαγ bγδ aβδ .
For thin shells, the deformation between S0 and S is fully characterized by aαβ and bαβ , and
their relation to Aαβ and Bαβ . Two important quantities that characterize the in-plane part of
the deformation are
I1 = Aαβ aαβ
(14)
7

and
I2 = det [Aαβ ] · det [aαβ ] .

(15)

They define the invariants of the√surface Cauchy Green tensors C = aαβ Aα ⊗ Aβ and B =
Aαβ aα ⊗ aβ . The quantity J = I2 determines the change of surface area between S0 and S.
Given the parametric derivative of aα , the so-called co-variant derivative of aα is given by
aα;β = (n ⊗ n) aα,β ,

(16)

i.e. as the change of aα along n. Introducing the Christoffel symbol Γγαβ := aγ · aα,β , we can
also write
aα;β = aα,β − Γγαβ aγ .
(17)
For general vectors, v = v α aα + v n, the parametric and co-variant derivative are considered
to agree, i.e. v ,α = v ;α . Therefore, v;α = v,α , n;α = n,α , (v α aα );β = (v α aα ),β and
α = v α + Γα v γ ,
v;β
,β
βγ

(18)

due to (10). Likewise definitions follow for the reference surface.
The variation and linearization of the above quantities can be found in Sauer and Duong (2015).

2.2

Quasi-static equilibrium

For negligible inertia, equilibrium of the thin shell is governed by the field equation (Steigmann,
1999)
T α;α + f = 0 ,
(19)
where f is a source term and
T α = N αβ aβ + S α n

(20)

is the traction vector acting on the surface normal to aα . It can be related to the stress tensor
σ = N αβ aα ⊗ aβ + S α aα ⊗ n

(21)

via Cauchy’s formula T α = σ T aα , or equivalently T = σ T ν with ν = να aα and T = να T α . In
Eqs. (20) and (21),
N αβ = σ αβ + bαγ M γβ ,
(22)
S α = −M βα;β
are the in-plane and shear stress components acting on the cross section. The stress and bending
moment components σ αβ and M αβ are given through constitution. Eqs. (19) and (22) are a
consequence of momentum balance (Sauer and Duong, 2015). At the boundary of the surface,
∂S, the boundary conditions
x = ϕ̄
on ∂x S ,
t

on ∂t S ,

= t̄

(23)

on ∂m S

mτ = m̄τ

can be prescribed. Here, mτ is the bending moment component parallel to boundary ∂S. For
Kirchhoff-Love shells, bending moments perpendicular to boundary ∂S, denoted mν , affect the
boundary traction. Therefore the effective traction
t := T − (mν n)0
is introduced (Sauer and Duong, 2015). We will consider mν = 0 in the following.
8

(24)

2.3

Constitution

The focus here is on the quasi-static behavior of lipid bilayers, which can be described in the
framework of hyperelasticity. For thin shells the stored energy function (per reference surface
area) takes the form W = W (aαβ , bαβ ), such that (Steigmann, 1999; Sauer and Duong, 2015)
2 ∂W
,
J ∂aαβ
1 ∂W
.
=
J ∂bαβ

σ αβ =
M αβ

(25)

From this we can then evaluate the in-plane stress N αβ and the shear S α according to (22). For
convenience, we further define τ αβ := Jσ αβ , M0αβ := JM αβ and N0αβ := JN αβ . The bending
behavior of lipid bilayers is commonly described by the bending model of Helfrich (1973)
w = k (H − H0 )2 + k ∗ κ ,

(26)

which is an energy density per current surface area. Here k is the bending modulus, k ∗ is the
Gaussian modulus and H0 denotes the so-called spontaneous curvature caused by the presence
of certain proteins. Based on (26), we consider the following two constitutive models:
2.3.1

Area-compressible lipid bilayer

Combining the Helfrich energy with an energy resulting from the surface area change, we write
W =Jw+

K
(J − 1)2 .
2

(27)

Here a simple quadratic term for the compressible part is considered, since the area change of
lipid bilayers is very small before rupture occurs (typically |J − 1| < 4%). According to (25)
and (22), the stress and moment components then become

σ αβ = K (J − 1) + k ∆H 2 − k ∗ κ aαβ − 2 k ∆H bαβ ,

M αβ = k ∆H + 2 k ∗ H aαβ − k ∗ bαβ ,
(28)

N αβ = K(J − 1) + k ∆H 2 aαβ − k ∆H bαβ ,
where ∆H := H − H0 .
Remark: Here, k and k ∗ are material constants that have the units strain energy per current
surface area. In principle, also Jk and Jk ∗ could be regarded as material constants that now
have the units strain energy per reference surface area. This alternative would lead to different
expressions for σ αβ and N αβ .
2.3.2

Area-incompressible lipid bilayer

Since K is usually very large for lipid-bilayers, one may as well consider the surface to be fully
area-incompressible. Using the Lagrange multiplier approach, we now have
W = Jw+qg ,

(29)

g := J − 1 = 0

(30)

where the incompressibility constraint

9

is enforced by the Lagrange multiplier q, which is an independent variable that needs to be accounted for. Physically, it corresponds to a surface tension. The stress and moment components
now become

σ αβ = q + k ∆H 2 − k ∗ κ aαβ − 2 k ∆H bαβ ,

M αβ = k ∆H + 2 k ∗ H aαβ − k ∗ bαβ ,
(31)

N αβ = q + k ∆H 2 aαβ − k ∆H bαβ ,
which is identical to (28) for q = Kg.
As K becomes larger and larger both models approach the same solution. So from a physical
point of view it may not make a big difference which model is used. Computationally, model
(27) is easier to handle but can become inaccurate for large K as is shown in Sec. 5. In analytical
approaches, often (29) is preferred as it usually simplifies the solution. Examples for the latter
case are found in Baesu et al. (2004) and Agrawal and Steigmann (2009); the former case is
considered in the original work of Helfrich (1973).
2.3.3

Model properties

In both models, the membrane part only provides bulk stiffness, but lacks shear stiffness. For
quasi-static computations the model can thus become unstable and should be stabilized, as is
discussed in Sec. 3. Interestingly, the bending part of the Helfrich model can contribute an
in-plane shear stiffness, which is shown in the following.
We first introduce the surface tension γ of the surface as the average trace of the stress tensor,
giving
γ := 21 σ : i = 12 Nαα .
(32)
For both (28) and (31) we find
γ = q − k H0 ∆H ,

(33)

where q = Kg in the former case. It can be seen that for H0 6= 0, the bending part contributes
to the surface tension. The surface tension is therefore not given by the membrane part alone
(Rangamani et al., 2014). For the compressible case, the effective bulk modulus can then be
determined from
∂γ
(34)
Keff :=
,
∂J
i.e. as the change of γ w.r.t J. We find
Keff = K + k H0 H/J ,

(35)

since ∂H/∂J = −H/J. Likewise we can define the effective shear modulus from
µeff := J aαγ

αβ
∂Ndev
aβδ ,
2 ∂aγδ

(36)

i.e. as the change of the deviatoric stress w.r.t to the deviatoric deformation (characterized by
aγδ /J). The deviatoric in-plane stress is given by
αβ
Ndev
:= N αβ − γ aαβ .

(37)

We find
αβ
Ndev
= k ∆H H aαβ − bαβ

10



(38)

for both (28) and (31). Evaluating (36) thus gives

µeff = Jk 3H 2 − 2HH0 − κ /2 .

(39)

The model therefore provides stabilizing shear stiffness if 3H 2 > 2HH0 + κ. Since this is not
always the case (e.g. for flat surface regions), additional shear stabilization should be provided
in general. This is discussed in Sec. 3. The value of µeff is discussed in detail in the examples
of Sec. 5. It is shown that µeff can sufficiently stabilize the problem such that no additional
shear stabilization is needed. It is also shown that µeff does not necessarily need to be positive
to avoid instabilities. Geometric stiffening, arising in large deformations, can also stabilize the
shell.

2.4

Weak form

The computational solution technique proposed here is based on the weak form governing the
mechanics of the lipid bilayer. Consider a kinematically admissible variation of the surface
denoted by δx ∈ V. Such a variation of S causes variations of aα , n, aαβ and bαβ . Contracting
field equation (19) with δx, integrating over S and applying Stokes’ theorem leads to the weak
form (Sauer and Duong, 2015)
G = Gint − Gext = 0 ∀ δx ∈ V ,
with

(40)

Z

Gint
Gext

Z
1
αβ
=
δaαβ σ da +
δbαβ M αβ da ,
S 2
S
Z
Z
Z
=
δx · f da +
δx · t ds +
δn · mτ ν ds .
S

∂t S

(41)

∂m S

Denoting the in-plane and out-of-plane components of δx by wα and w, such that δx := wα aα +
w n, we find that
δaαβ = wα;β + wβ;α − 2w bαβ .
(42)
Thus, the first part of Gint can be split into in-plane and out-of-plane contributions as
Z
1
out
δaαβ σ αβ da = Gin
σ + Gσ ,
S 2
with
Gin
σ

Z
=

and
Gout
σ =−

(43)

wα;β σ αβ da

(44)

w bαβ σ αβ da .

(45)

S

Z
S

In principle – although not needed here – the second part of Gint can also be split into in-plane
and out-of-plane contributions (Sauer and Duong, 2015). In the area-incompressible case, we
additionally have to satisfy the weak form of constraint (30),
Z
Gg =
δq g dA = 0 ∀ δq ∈ Q ,
(46)
S0

where Q is an admissible space for the variation of Lagrange multiplier q.

11

3

Liquid shell stabilization

As noted above, the system is unstable for quasi-static computations. There are two principal ways to stabilize the system without modifying/affecting the original problem. They are
discussed in the following two sections and then summarized in Sec. 3.3.

3.1

Adding stiffness

αβ
Firstly, the system can be stabilized by adding a stabilization stress σsta
to σ αβ in order to
provide additional stiffness. This stress can be defined from a (convex) shear energy or from
numerical viscosity. An elegant and accurate way to stabilize the system is to add the stabilization stress only to the in-plane contribution (44) while leaving the out-of-plane contribution
(45) unchanged. The advantage of this approach is that the out-of-plane part, responsible for
the shape of the bilayer, is not affected by the stabilization, at least not in the continuum limit
of the surface discretization. There are several different ways to define the stabilization stress,
which we will group into two categories. An overview of all the options is then summarized in
Tab. 1.

3.1.1

In-plane shear and bulk stabilization

The first category goes back to Sauer (2014), who used it to stabilize liquid membranes governed
by constant surface tension. The stabilization stress for such membranes requires shear and bulk
contributions. Those are given for example by the stabilization stress

αβ
(47)
σsta
= µ/J Aαβ − aαβ ,
based on numerical stiffness, and

αβ
αβ ,
σsta
= µ/J aαβ
pre − a

(48)

αβ at the preceding computational
based on numerical viscosity. Here aαβ
pre denotes the value of a
step. These stabilization stresses are then only included within Eq. (44) and not in Eq. (45),
and the resulting two stabilization schemes are denoted ‘A’ (for (47)) and ‘a’ (for (48)) following
Sauer (2014). This reference shows that scheme ‘a’ is highly accurate and performs much better
than scheme ‘A’. It is also shown that applying the stabilization stresses (47) and (48) only to
the in-plane part is much more accurate than applying it throughout the system (i.e. in both
Eqs. (44) and (45)), which we denote as schemes ‘A-t’ and ‘a-t’.

3.1.2

Sole in-plane shear stabilization

If the surface tension is not constant, as in the lipid bilayer models introduced above, only
shear stabilization is required. Therefore, the following new category of stabilization schemes
is defined. Consider a split of the surface Cauchy-Green tensor into dilatational and deviatoric
b where C
b := J −1 C describes only the deviatoric deformation (since
parts, such that C = J C,
b = 1). The stored energy function of an elastic membrane can then be defined for example
det C
by
 µ

K 2
(49)
W =
J − 1 − 2 ln J +
Ib1 − 2 ,
4
2

12

b The first term in (49) captures purely dilatoric
where Ib1 = I1 /J is the first invariant of C.
deformations, while the second part captures purely deviatoric deformations. The formulation
is analogous to the 3D case described for example in Wriggers (2008). Since the bilayer energy
introduced in Sec. 2.3 already contains a bulk part, we now only need to consider the contribution
Wsta =


µ b
I1 − 2
2

(50)

to derive the new stabilization scheme. From (25) we then find the stabilization stress
µ  αβ I1 αβ 
A − a
,
J2
2

αβ
σsta
=

(51)

or

µ  αβ I1 αβ 
(52)
A − a
.
J
2
As before, this stress will only be applied to Eq. (44) and not in Eq. (45), even though it has
been derived from a potential and should theoretically apply to both terms. Following earlier
nomenclature we denote this scheme by ‘A-s’. Replacing Aαβ by aαβ
pre in (52) gives
αβ
τsta
=

αβ
τsta
=

with J ∗ :=

I1∗ αβ 
µ  αβ
a
,
−
a
J ∗ pre
2

(53)

q

αβ
∗
det aαβ det apre
αβ and I1 := apre aαβ , which is an alternative shear-stabilization

scheme based on numerical viscosity. We denote it ‘a-s’. If stresses (52) and (53) are applied
throughout the system (i.e. to both (44) and (45)), we denote the corresponding schemes ‘A-st’
and ‘a-st’.
If the shell is (nearly) area-incompressible the two stabilization methods of Sec. 3.1.1 and 3.1.2
can behave identical, as can be seen by example 5.1.
3.1.3

Relation to physical viscosity

Schemes (48) and (53) are related to physical viscosity. Considering (near) area-incompressibility
(J ∗ = J = 1), the viscous stress for a Newtonian fluid is given by
αβ
σvisc
= −ν ȧαβ

(54)

(Aris, 1989; Rangamani et al., 2013, 2014), where ȧαβ = −aαγ ȧγδ aδβ . Considering the first
order rate approximation

1  αβ
(55)
ȧαβ ≈
a − aαβ
pre ,
∆t
and a small time step (I1∗ ≈ 2), immediately leads to expressions (48) and (53) with
ν = µ ∆t .

3.2

(56)

Normal projection

The second principal way to stabilize the system consists of a simple projection of the formulation onto the solution space defined by the normal surface direction. We apply this step directly
to the discretized formulation as was proposed by Sauer (2014). According to this, for the given

13

discrete system of linear equations for displacement increment ∆u given by K ∆u = −f , the
reduced system for increment ∆ured = P ∆u is simply obtained as
Kred ∆ured = −fred , Kred := P K PT , fred := P f ,
where



nT1

0T

···

0T






P := 



0T
..
.

nT2
..
.

···
..
.

0T
..
.








0T

0T

···

nTnno

(57)

(58)

is a projection matrix defined by the nodal normal vectors nI . Since this method can lead
to distorted FE meshes, a mesh update can be performed by applying any of the stabilization
techniques discussed above. If this is followed by a projection step at the same load level, a
dependency on parameter µ is avoided.
For NURBS discretizations, the computation of the normal nI corresponding to control point I
is not trivial, since the control points generally do not lie on the surface. The control points can
be projected onto the actual surface in order to evaluate the normal at the projection point. As
a simplification of this approach one can also work with the approximate normal obtained at the
current location of the initial projection point. The difference between the two approaches is
very small as Fig. 1 shows. There are minor differences in the locations of the projection points

a.

b.

Figure 1: Search direction nI for stabilization scheme ‘P’: a. Considering the normals at the
projected control points; b. Considering the normals at the current location of the initially
projected control points.

where nI is evaluated. But nI itself does not change much. In the example considered here and
in Sauer (2014), no siginificant difference was therefore found in the numerical results. In order
to solve the reduced problem (57) appropriate boundary conditions are needed. Boundary
conditions can be applied before reduction or they can be suitably adapted to the reduced
system. This is also shown in Fig. 1. In this example (see Sec. 5.3) the vertical displacement at
the outer two and inner two rings of nodes is prescribed. After a mesh update at a given solution
step no further vertical displacement can therefore occur on those nodes. But the nodes can
still move radially. Therefore one can simply replace nI at those nodes by radial unit vectors.
Remark: For liquid droplets, stabilization scheme ‘P’ shows superior accuracy over all other
schemes (Sauer, 2014). However this is not the case for area-incompressible surfaces: For a
curved surface, non-zero values of ∆ured lead to a change in surface area. In regions of non-zero
14

curvature, ∆ured will therefore tend to be zero, yielding scheme ‘P’ ineffective. This is confirmed
by the results in Sec. 5.3.

3.3

Summary of the stabilization schemes

The nine stabilization schemes presented above are summarized in Tab. 1. They can be grouped
αβ
application of σsta

dependence

only in (44)

only on µ

both in (44) & (45)

only on µ

only in (44)

only on µ

both in (44) & (45)

only on µ

only in (44)

on µ and nt

both in (44) & (45)

on µ and nt

only in (44)

on µ and nt

a-st

αβ
stab. stress σsta
/µ

αβ
αβ
A − a /J

Aαβ − aαβ /J

Aαβ − 12 I1 aαβ /J 2

Aαβ − 12 I1 aαβ /J 2

αβ /J
aαβ
pre − a

αβ /J
aαβ
pre − a
 ∗2
1 ∗ αβ
aαβ
/J
pre − 2 I1 a
 ∗2
αβ
1 ∗ αβ
apre − 2 I1 a /J

both in (44) & (45)

on µ and nt

P

0

–

on nodal nI

class

scheme

A

A
A-t
A-s
A-st

a

a
a-t
a-s

P

Table 1: Summary of the stabilization schemes presented in Sec. 3.1.1 and 3.1.2
into three classes: A, a and P. The schemes of class A depend only on µ but require this value
to be quite low. The schemes of class a also depend on the number of computational steps, nt . If
this number is high the schemes provide stiffness without adding much stress. The shell is then
stabilized without modifying the solution much, even when µ is high. Scheme ‘P’ depends on
the nodal projection vector nI , which is usually taken as the surface normal. The performance
of the different stabilization schemes is investigated in the examples of Sec. 5.

4

FE formulation

This section presents a general finite element formulation for lipid bilayers using both material
models of Sec. 2.3. The formulation applies to general 3D surfaces. It is based on the solid shell
formulation of Duong et al. (2015).

4.1

FE approximation

The geometry of the reference surface and the current surface are approximated by the finite
element interpolations
X ≈ N Xe
(59)
and
x ≈ N xe ,

(60)

where N := [N1 1, ..., Nne 1] is a (3 × 3ne ) array containing the ne nodal shape functions
T T
NI = NI (ξ 1 , ξ 2 ) of element Ωe defined in parameter space. Xe := [X T
and
1 , ...., X ne ]

15

T T
e
1
xe := [xT
1 , ...., xne ] contain the ne nodal position vectors of Ω . In order to ensure C continuity of the surface, NURBS-based shape functions are used (Borden et al., 2011). The
tangent vectors of the surface are thus approximated by

∂X
≈ N,α Xe
∂ξ α

Aα =

(61)

and

∂x
≈ N,α xe ,
(62)
∂ξ α
From these, the normal vectors of the surface are then determined from (2) and (3). The
variation of x and aα are approximated in the same manner, i.e.
aα =

δx ≈ N δXe

(63)

δaα ≈ N,α δxe .

(64)

and
Based on these expressions, all the kinematical quantities discussed in Sec. 2.1 as well as their
variation (Sauer and Duong, 2015) can be evaluated.

4.2

Discretized weak form

Inserting the above interpolations into the weak form of Sec. 2.4 leads to the approximated
weak form
nel
nel
X
X
e
G≈
G =
(Geint − Geext ) .
(65)
e=1

e=1

The internal virtual work of each finite element is given by (Duong et al., 2015)

e
e
Geint = δxT
e fintσ + fintM ,

(66)

with the FE force vectors due to the membrane stress σ αβ and the bending moment M αβ
Z
e
fintσ
:=
σ αβ NT
(67)
,α aβ da ,
Ωe

and
e
fintM

Z
:=
Ωe


γ
T
M αβ NT
,αβ − Γαβ N,γ n da .

(68)

e
Following decomposition (43), fintσ
can be split into the in-plane and out-of-plane contributions
(Sauer et al., 2014)
e
e
e
finti
:= fintσ
− finto
,
Z
(69)
e
finto
:= −
σ αβ bαβ NT n da .
Ωe

The external virtual work for each FE due to the external loads f , t and mτ is given by (Duong
et al., 2015)

e
e
e
Geext = δxT
(70)
e fextf + fextt + fextm ,
with
e
fextf

Z
Ωe

e
fextt
:=

Z
∂t

e
fextm

NT f da ,

:=

NT t ds ,
Ωe

Z
:=
∂m Ωe

α
NT
,α ν mτ n ds .

16

(71)

For the complete linearization of these terms, see Sauer et al. (2014), Duong et al. (2015) and
Appendix A. The linearization requires the material tangents of σ αβ and M αβ . For the material
models in (28) and (31) these are given in Sauer and Duong (2015). For the stabilization schemes
in Sec. 3.1.1 the tangent is given in Sauer (2014), while the tangent for schemes ‘A-s’ and ‘A-st’
(see (52)) in Sec. 3.1.2 is


∂τ αβ
µ I1 αβ γδ
αβγδ
αβ γδ
αβ γδ
αβγδ
(72)
c
:= 2
.
=
a a − I1 a
−a A −A a
∂aγδ
J 2
For schemes ‘a-s’ and ‘a-st’, the terms I1 and J are simply substituted by I1∗ and J ∗ , respectively.
Here, aαβγδ is given in Sauer and Duong (2015).

4.3

Area constraint

For the area-incompressible case, a further step is to discretize the constraint and the corresponding Lagrange multiplier. For the latter, we write
q ≈ L qe ,

(73)

as in Eq. (60). Here L := [L1 , ..., Lme ] is a (1 × me ) array containing the me nodal shape
functions LI = LI (ξ 1 , ξ 2 ) of surface element Ωe , and qe := [q1 , ...., qme ]T contains the me nodal
Lagrange multipliers of the element. It follows that
δq ≈ L δqe ,

(74)

such that weak constraint (46) becomes
nel
X

Geg ,

(75)

e
Geg = δq T
e g ,

(76)

Gg ≈

e=1

where
with
ge

Z

LT g dA .

:=
Ωe0

(77)

The linearization of ge , needed for the following solution procedure, is provided in Appendix B.

4.4

Solution procedure

e , fe
e
e
e
e
The elemental vectors fintσ
intM , fextf , fextt , fextm and g are assembled into the global vectors
f and g by adding corresponding entries. The discretized weak form then reads

δxT f (x, q) + δqT g(x) = 0 ,

∀ δx ∈ V h & δq ∈ Qh ,

(78)

where x, q, δx and δq are global vectors containing all nodal deformations, Lagrange multipliers
and their variations. V h and Qh are corresponding discrete spaces. Eq. (78) is satisfied if f = 0
and g = 0 at nodes where no Dirichlet BC apply. These two nonlinear equations are then solved
with Newton’s method for the unknowns x and q. We note that the discretization of x and
q should satisfy the LBB-condition (Babuška, 1973; Bathe, 1996). For that, we consider here
C 1 -continuous, bi-quadratic NURBS interpolation for x and C 0 -continuous, bi-linear Lagrange
interpolation for q. If no constraint is present (like in model (27)), the parts containing q and
g are simply skipped. A comparison between the different models ((27) and (29)) is presented
in Sec. 5.
17

4.5

Rotational constraints

To constrain rotations, we add the constraint potential
Z

Πn =
(n − n̄) · (n − n̄) dS
L0 2

(79)

to the shell formulation. This approach can be used to apply rotations at boundaries, to fix
rotations at symmetry boundaries, and to equalize normals at patch boundaries. The particular
boundary under consideration is denoted as L0 in the reference configuration.  is a penalty
parameter. The variation, linearization and FE discretization of (79) is discussed in Duong
et al. (2015).

4.6

Normalization

For a numerical implementation, the above expressions need to be normalized. We normalize
the geometry and deformation by some length scale L, i.e. X̄ = X/L and x̄ = x/L, where the
bar indicates non-dimensional quantities, and non-dimensionalize all the kinematics based on
this. The material parameters are chosen to be normalized by parameter k, which has the unit
[force × length]. The non-dimensional material parameters thus are
k̄ = 1 ,
k̄ ∗ = k ∗ /k ,
K̄ = K L2 /k ,
µ̄ =

µ L2 /k

(80)

,

¯ =  L/k .
With the chosen normalization parameters k and L, the normalization of stress and moment
components become3
q̄
= q L2 /k ,
σ̄ αβ = σ αβ L2 /k ,
M̄ αβ

=

M αβ

(81)

L/k ,

while the normalization of the loading follows as
f̄

= f L3 /k ,

t̄

= t L2 /k ,

(82)

m̄τ = mτ L/k .

5

Numerical examples

To illustrate the performance of the proposed finite element model, four examples – marked
by increasing computational complexity – are presented here. The first three examples have
analytical solutions that are used for model verification.

5.1

Pure bending and stretching of a flat strip

The first example considers the bending of a flat strip, with dimension S × L, by applying the
rotation Θ/2 at the ends of the strip as is shown in Fig. 2a. Further, a uniform stretch, with
3

supposing that the parameters ξ α carry units of length, so that aαβ and aαβ become dimensionless

18

a.

b.

Figure 2: Pure bending: a. initial FE configuration and boundary conditions (for S = πL
discretized with m = 8 elements); b. current FE configuration for an imposed rotation and
stretch of Θ/2 = π/3 and λ2 = 1.5. The color shows the relative error in the mean curvature,
defined as eH := 1 − HFE /Hexact , considering model (29) without any stabilization.

magnitude λ2 , is applied as shown in the figure. The remaining boundaries are supported as
shown. In particular, the rotation along the boundaries at Y = 0 and Y = L is not constrained.
We will study the problem by examining the support reactions M = M (Θ, λ2 ) (the distributed
moment along X = 0) and N = N (Θ, λ2 ) (the traction along Y = 0). The analytical solution
for this problem is given by Sauer and Duong (2015). According to this, the strip deforms into
a curved sheet with dimension s × ` = λ1 S × λ2 L and constant mean curvature
H=

κ1
,
2

(83)

where λ1 = s/S, λ2 = `/L and κ1 = Θ/s are the in-plane stretches and the out-of-plane
curvature of the strip. Further,
 −2

λ1
0
αβ
(84)
[a ] =
0 λ−2
2
and
[bαβ ]


=

κ1 λ−2
0
1
0
0


.

(85)

With this, the in-plane stress components become
N11 = q − k H 2 ,
N22 = q + k H 2

(86)

both for the area-incompressible model of (29) and the compressible model of (27) with q =
K(J − 1). For the considered boundary conditions, N11 = 0, so that
q = k H 2,

(87)

and we thus have the support reaction (per current length) N := N22 = 2k H 2 along Y = 0
and Y = `. Per reference length this becomes N0 = λ1 N . The bending moment necessary to
support the applied rotation (along X = 0 and X = πR) becomes (Sauer and Duong, 2015)
M = kH
19

(88)

per current length of the support (or M0 = λ2 M per reference length). If the special case
k ∗ = −k/2 (Canham, 1970) is considered, there is no bending in the Y -direction.
For the area-incompressible model of (29), we have λ1 = 1/λ2 . For the area-compressible case
according to model (27), we can determine λ1 from (87) with J = λ1 λ2 , giving
λ1 =

i
1hk 2
H +1 .
λ2 K

(89)

The two cases are solved numerically using the computational setup shown in Fig. 2. The FE
mesh is discretized by m elements along X. The parameter t is introduced to apply the rotation
Θ = tπ/6 and stretch λ2 = 1 + t/2 by increasing t linearly from 0 to 1 in nt steps, where nt
is chosen as a multiple of m. The mean curvature then follows as H = Θ/(2λ1 S). For the
unconstrained case, λ1 is then the solution of the cubic equation
λ31 λ2 − λ21 −

Θ2
=0,
4K̄

(90)

with K̄ = KS 2 /k. Numerically, the rotation is applied according to (79) considering the penalty
parameter  = 100 nx k/L. Fig. 3 shows the FE solution and analytical solution for M0 (t) and
N0 (t), normalizing M by k/L and N by k/L2 .
0.8

0.5
FE M0(t)

0.6

exact N0(t)

0.5

exact M0(t)
M0 [k/L] or N0 [k/L2]

0.7

FE N0(t)

0.4
0.3

0

0

M [k/L] or N [k/L2]

exact M0(t)

0.2

0.4

FE M0(t)
exact N0(t)
FE N0(t)

0.3

0.2

0.1
0.1

a.

0
0

0.2

0.4
0.6
parameter t

0.8

1

b.

0
0

0.2

0.4
0.6
parameter t

0.8

1

Figure 3: Pure bending: distributed boundary moment M0 (t) and normal traction N0 (t) as
obtained analytically and computationally for a. the area-incompressible model (29) and b. the
area-compressible model (27) for K̄ = 2.5.

Next, the accuracy of the different stability schemes is studied in detail by examining the L2 error of the solution, defined by
s
Z
1
(91)
L2 :=
kuexact − uFE k2 dA ,
SL S0
and the error in M and N , defined by
EM N :=

|Mexact − MFE | |Nexact − NFE |
+
,
Mexact
Nexact

(92)

where MFE and NFE are the computed mean values along the respective boundaries. The
first error is a measure of the kinematic accuracy, while the second is a measure of the kinetic
20

0

0

10

10

−1

−1

10
L2 error of U

L2 error of U

10

−2

10

stab. A−s
stab. A−t
stab. A
stab. A−st
stab. a−s
stab. a−t
stab. a
stab. a−st

−3

10

−4

10

−5

10

0

a.

−5

10
Number of elements

10

2

10

−10

10
stabilization parameter µ [k/L2]

b.

2

0

10

2

10

1

1

10

10

stab. A−s
stab. A−t
stab. A
stab. A−st

0

0

10

EMN

EMN

10
stab. A−s
stab. A−t
stab. A
stab. A−st
stab. a−s
stab. a−t
stab. a
stab. a−st

−1

10

−2

10

−3

c.

−3

10

−4

10

10

−2

10

10

1

10

stab. A−s
stab. A−t
stab. A
stab. A−st

0

10

−1

10

−2

10

−3

1

10
Number of elements

10

2

10

d.

−10

10
stabilization parameter µ [k/L2]

0

10

Figure 4: Pure bending: accuracy for the area-incompressibile model (29): a. L2 -error vs. m
considering stabilization classes A and a with µ̄ = 10 and nt = 12.5 m; b. L2 -error vs. µ
considering stabilization class A with m = 32; c.–d. same as a.–b., but now for error EM N .
Considered is Θ/2 = π/3 and λ2 = 1.5.

accuracy. Fig. 4 shows the two errors for the area-incompressible model of Eq. (29). Looking at
the L2 -error, schemes ‘A-t’, ‘A-st’, ‘a-t’ and ‘a-st’ perform best. In case of error EN M , schemes
‘a’ and ‘a-s’ perform best. Class A generally converges with µ, but it may not converge with the
number of elements for high values of µ. Interestingly, the L2 -error of scheme ‘A-t’ and ‘A-st’
is not affected by µ, as schemes ‘A’ and ‘A-s’ are. For sufficiently low µ (for m = 32 about
µ̄ < 10−3 ), the accuracy of class A (both in L2 and EM N ) reaches that of class a and then only
improves with mesh refinement. Class A with low µ may even surpass class a with high µ. But
generally, class a is more accurate and robust (as µ does not need to be very small). There is
no clear favourite in class a for this test case.
Fig. 5 shows the two errors for the area-compressible model of Eq. (27) considering K̄ = 2.5.
In case of the L2 -error, scheme ‘a’ performs best, while for error EM N , scheme ‘a-s’ is best. As
before, class A is poor for large µ, but reaches (and may surpass) the accuracy of class a at
some µ depending on m.
As the plots show, not a single stabilization scheme stands out here and the accuracy depends
both on the model and the error measure. In general, all schemes are suitable to solve the
problem. If class A is used, the value of µ needs to be suitably low. For class a even large
values for µ can be used. In this example it is even possible to set µ = 0 in the code. This works
21

0

0

10

10

−1

−1

10
L2 error of U

L2 error of U

10

−2

10

stab. A−s
stab. A−t
stab. A
stab. A−st
stab. a−s
stab. a−t
stab. a
stab. a−st

−3

10

−4

10

−5

10

0

a.

−5

10
Number of elements

10

2

10

−10

10
stabilization parameter µ [k/L2]

b.

2

0

10

2

10

1

1

10

10

stab. A−s
stab. A−t
stab. A
stab. A−st

0

0

10

EMN

EMN

10
stab. A−s
stab. A−t
stab. A
stab. A−st
stab. a−s
stab. a−t
stab. a
stab. a−st

−1

10

−2

10

−3

c.

−3

10

−4

10

10

−2

10

10

1

10

stab. A−s
stab. A−t
stab. A
stab. A−st

0

10

−1

10

−2

10

−3

1

10
Number of elements

10

2

10

d.

−10

10
stabilization parameter µ [k/L2]

0

10

Figure 5: Pure bending: accuracy of the area-compressible model (27): a. L2 -error vs. m
considering stabilization classes A and a with µ̄ = 10 and nt = 12.5 m; b. L2 -error vs. µ
considering stabilization class A with m = 32; c.–d. same as a.–b., but now for error EM N .
Considered is K̄ = 2.5, Θ/2 = π/3 and λ2 = 1.5.

since the effective shear stiffness according to (39) is positive here, i.e. µeff = 3JkH 2 /2 > 0. For
other problems µeff can be negative, and stabilization is required.

5.2

Inflation of a sphere

The second example considers the inflation of a spherical cell. Contrary to the previous example,
the FE mesh now also contains interfaces between NURBS patches. Since the surface area
increases during inflation, potential (27) is considered. For this model, the in-plane traction
component, given in (28), is
N αβ = Na aαβ + Nb bαβ ,
(93)
with

Na := k ∆H 2 + K (J − 1) ,
Nb := −k ∆H .

(94)

The initial radius of the sphere is denoted by R, the initial volume is denoted by V0 = 4πR3 /3.
The cell remains spherical during inflation, so that we can obtain an analytical reference solution.

22

The current radius during inflation shall be denoted by r, the current volume by V = 4πr3 /3.
Considering the surface parameterization


r cos φ sin θ
x(φ, θ) =  r sin φ sin θ  ,
(95)
−r cos θ
we find
[aαβ ]

1
= 2
r

"

1/ sin2 θ 0
0

1

#
,

(96)

bαβ = −aαβ /r and H = −1/r for this example. The traction vector T = να T α on a cut ⊥ ν
thus becomes
T = (Na − Nb /r) ν + S α να n
(97)
according to (20). The in-plane component Tν := Na − Nb /r must equilibrate the current
pressure according to the well-known relation
p=

2Tν
.
r

(98)

We can thus establish the analytical pressure-volume relation
 1

1
1
2
p̄(V̄ ) = 2H̄0 V̄ − 3 − 2H̄02 V̄ − 3 + 2K̄ V̄ 3 − V̄ − 3 ,

(99)

normalized according to the definitions p̄ := pR3 /k, V̄ := V /V0 , H̄0 := H0 R and K̄ := KR2 /k.
Fig. 6 shows the computational setup of the problem. The computational domain consists of

a.

b.

Figure 6: Sphere inflation: a. initial FE configuration and boundary conditions (for mesh
m = 8); b. current FE configuration for an imposed volume of V̄ = 2 compared to the initial
configuration; the colors show the relative error in the surface tension Tν .
a quarter sphere discretized with four NURBS patches. The quarter sphere contains 3m2 /2
elements where m is the number of elements along the equator of the quarter sphere. At the
boundaries and at the patch interfaces C 1 -continuity is enforced using (79) with  = 4000mk/R.
The area bulk modulus is taken as K = 5k/R2 , while k ∗ is taken as zero. Two cases are
considered: H0 = 0 and H0 = 1/R. Fig. 7 shows that the computational p(V )-data converge to
the exact analytical result of (99). Here the pressure error
ep =

|pexact − pFE |
pexact
23

(100)

−2

8

10

−3

10

4
error in p

internal pressure [ k/R3 ]

6

2
0
exact H0 = 0

−2

−5

FE H0 = 0

10

FE H0 = 1

a.

A−t
a−t
all other schemes

exact H0 = 1

−4
−6
0.5

−4

10

1
1.5
enclosed volume [ V/V0 ]

−6

10

2

0

10

b.

1

10
mesh parameter m

2

10

Figure 7: Sphere inflation: a. pressure-volume relation; b. FE convergence for the different
stabilization schemes.

is examined for H0 = 1/R and V̄ = 2 considering the 9 stabilization schemes of Tab. 1 with
µ̄ = 0.01 for class A and µ̄ = 1 and nt = 5m for class a. For schemes ‘A’, ‘A-s’, ‘A-st’, ‘a’, ‘a-s’,
‘a-st’ and ‘P’ this error converges nicely (and is indistinguishable in the figure). Only schemes
‘A-t’ and ‘a-t’ behave significantly different. They introduce further errors that only converge if
µ is decreased or nt is increased. The reason why all other schemes have equal error, is that here
the error is actually determined by the penalty parameter  used within patch constraint (79).
The error stemming from the stabilization methods (apart from ‘A-t’ and ‘a-t’) is insignificant
compared to that. It is interesting to note that ‘A-st’ and ‘a-st’ perform much better than
‘A-t’ and ‘a-t’, even though no shear is present in the analytical solution. ‘A-st’ and ‘a-st’ can
therefore be considered as the best choices here, since they are the most efficient schemes to
implement.
We finally note that for a sphere µeff = JkH(H − H0 ), where H = −1/r. Thus µeff > 0 for
H < H0 , which is the case here.

5.3

Drawing of a tube

Tube-like shapes are one of the most common non-trivial shapes in biological cell membranes.
They can be observed in the endoplasmic reticulum (Shibata et al., 2006, 2009; Hu et al., 2009)
and mitochondria (Fawcett, 1981; Griparic and van der Bliek, 2001; Shibata et al., 2009), where
individual tubules or a complex network of tubules co-exist. The tubes can also be formed when
a membrane is pulled by actin or microtubule polymerization, or by molecular motors attached
to membranes but walking on cytoskeletal filaments (Terasaki et al., 1986; Waterman-Storer
and Salmon, 1998; Koster et al., 2003; Itoh et al., 2005). These situations can be modeled by
means of a lateral force that acts on a membrane. An analytical solution for this problem has
been obtained by Derényi et al. (2002) by assuming axisymmetry and infinitely long tubes,
which we use to verify our finite element formulation. The dynamics of tether formation and
relaxation have been also studied by Rahimi and Arroyo (2012).
To simulate the tube drawing process we consider the setup shown in Fig. 8. The cell membrane
is modeled as a circular, initially flat disc with initial radius L. The effect of the surrounding
membrane is captured by the boundary tension σ (w.r.t. the current boundary length). The
surface is described by material model (27). L and k are used for normalization. The remaining

24

a.

b.

Figure 8: Tube drawing: a. boundary conditions and b. coarse FE mesh of the initial configuration. The dots show the control points of the NURBS discretization.

material parameters are chosen as k̄ = −0.7 k and K = 20,000 k/L2 . H0 is taken as zero.
We consider σ ∈ {100, 200, 400, 800}k/L2 . Stabilization scheme ‘A-s’ is considered with µ =
0.1 k/L2 . The shell is clamped at the boundary, but free to move in the in-plane direction.
The traction t = σν is imposed and applied numerically via (71.2). Even though t is constant
e
changes and has to be linearized.
during deformation, the boundary length ds appearing in fextt
This is, for example, discussed in Sauer (2014). At the center, the displacement u is imposed
on the initially flat, circular surface. To prevent rigid rotations, selected nodes are supported
as shown in Fig. 8b.
Fig. 8b also shows one of the chosen finite element discretizations of the initial configuration.
Quadratic, NURBS-based, C 1 -continuous finite elements are used. For those, the outer ring of
control points lies outside of the mesh; their distance to the center is therefore slightly larger
than L. A finer discretization is chosen at the center, where the tube is going to form. The
chosen NURBS discretization degenerates at the center, such that the C 1 -continuity is lost
there. It is regained if displacement u is applied not only to the central control point but also to
the first ring of control points around the center. This ensures that the tangent plane remains
horizontal at the tip. Likewise, a horizontal tangent is enforced at the outer boundary by fixing
the height of the outer two rings of control points.
Fig. 9 shows the deformed surface at u = L, considering different values of the far-field surface
tension σ. The surface tension affects the slenderness of the appearing tube. Derényi et al.
(2002) showed from theoretical considerations4 that the tube radius is
r
1 k
(101)
a=
,
2 σ
while the steady force during tube drawing is
√
P0 = 2π σ k .
These values are confirmed by our computations, which is illustrated in Fig. 10.

4

Assuming that the tube is sufficiently long and can be idealized by a perfect cylinder.

25

(102)

Figure 9: Tube drawing: Results for σ ∈ {100, 200, 400, 800} k/L2 (bottom to top); colorscale
shows the mean curvature H normalized by L−1 .

26

0

1.2

10

µ = 0.1
µ = 0.4/m
µ=0
analytical

1
−1

0.8

error in P

drawing force P/P0

10

0.6

−2

10

0.4

0
0

−3

σ = 800
σ = 400
σ = 200
σ = 100

0.2

10

a.

20
30
40
prescribed displacement u/a

50

10

−4

10

60

0

10

b.

1

2

10
10
mesh parameter m

3

10

Figure 10: Tube drawing: a. load-displacement curve; b. FE convergence.

The left inset shows the force-displacement relation during drawing. Oscillations appear in the
numerical solution due to the mesh discretization error. They are more pronounced for more
slender tubes, as the black curve in Fig. 10a shows. They disappear upon mesh refinement, as
the solution converges. The convergence of P0 for σ = 200k/L2 and u = L (= 28.28a) is shown
in Fig. 10b by examining the error
e(P0FE ) :=

|P0ref − P0FE |
,
P0ref

(103)

where P0ref is the FE solution for m = 256 and µ = 0. The mesh sequence given in Tab. 2 is
used for the convergence study. For the convergence study, the following stabilization cases are
m
4
8
16
32
64
128
256

nel = 4 m2
64
256
1024
4096
16384
65336
262144

nno = 4m(m + 1)
80
288
1088
4224
16640
66048
263168

µ̄
1/10
1/20
1/40
1/80
1/160
1/320
1/640

Table 2: tube drawing: computational parameters (number of surface elements nel , number of
control points nno and stability parameter µ̄) for the convergence study of Fig. 10b.
considered:
1. scheme ‘A-s’ with fixed µ̄ = 0.1,
2. scheme ‘A-s’ with varying µ̄, as specified in Tab. 2,
3. scheme ‘P’ using the solution of case 1 as initial guess,
4. scheme ‘P’ using the solution of case 2 as initial guess,
5. no stabilization (µ = 0), using the solution of case 2 as initial guess.
It turns out that scheme ‘P’ does not at all improve the results obtained by scheme ‘A-s’,
probably due to the issue noted in the remark of Sec. 3.2. Case 1 does not converge below an
27

error of about 0.6%, which reflects the error caused by µ̄ = 0.1. Case 5 works due to the inherent
shear stiffness of the Helfrich model given in (39). For the cylindrical part, µ̄eff = 3H 2 /2 > 0.
With H = −1/(2a) follows µ̄eff = 3/(2a2 ). At the tip, the principal curvatures are equal
(κ1 = κ2 ), so that µ̄eff = κ21 > 0. At the tip the curvature is almost twice as large as the
cylinder curvature (i.e. κ1 ≈ −2/a), so that µ̄eff ≈ 4k/a2 . In the cap, µ̄eff ranges in between
these two extreme values, and so µeff is always positive. The initial flat disk has µeff = 0, so
that stabilization is needed for the initial step. In all cases, the error reported in Fig. 10b is
assessed by comparison to the finest mesh of case 2. From this we can find that the analytical
solution itself has an error of about 0.2%, due to its underlying assumptions.

5.4

Cell budding

The last example considers the budding of spherical cells. The example is used to demonstrate
the capabilities of the proposed computational model, in particular in the context of non-trivial
and non-axisymmetric deformations.
5.4.1

Motivation

It is known that protein adsorption can mediate shape changes in biological membranes (Zimmerberg and Kozlov, 2006; McMahon and Gallop, 2005; Kozlov et al., 2014; Shi and Baumgart,
2015). The lipid membrane deforms whenever its curvature is incompatible with the inherent
structure of a protein, giving rise to a spontaneous curvature. Such protein-induced spontaneous curvature is common in many important biological phenomena such as endocytosis
(Buser and Drubin, 2013; Kukulski et al., 2012; Peter et al., 2004), cell motility (Keren, 2011)
and vesicle formation (Gruenberg and Stenmark, 2004). For example, in endocytosis, a clathrin
protein coat is formed on the membrane which is incompatible with the membrane curvature,
thus driving the out-of-plane shape changes, finally leading to fission. Moreover, another set of
curvature generating proteins, the so-called BAR proteins (Peter et al., 2004; Kukulski et al.,
2012), play an important role in modulating shape changes leading to fission in the later stages
of endocytosis.
In this example, we use the proposed finite element formulation to study the shape changes in
membranes due to protein-induced spontaneous curvature. This spontaneous curvature usually
depends on the concentration of proteins adsorbed onto the membrane. In obtaining the shapes,
we assume that any proteins that are bound to the membrane do not diffuse and are concentrated
in a specific region (Walther et al., 2006; Agrawal and Steigmann, 2009; Karotki et al., 2011;
Kishimoto et al., 2011; Rangamani et al., 2014). Such processes are common in all endocytosis
related phenomena, where it is known that the clathrin protein coat does not diffuse in the
membrane once adsorbed. Therefore, our example aims at helping to understand the shapes
arising in the early stages of endocytosis.
5.4.2

Computational setup

For our example, we consider a hemi-spherical cell with initial radius R; the cell surface is
clamped at the boundary, but free to expand radially as is shown in Fig. 11. On the tip
of the cell, within the circular region of radius 0.2R, a constant spontaneous curvature H̄0
is prescribed. Unless otherwise specified, model (27) is used with the material parameters
k̄ ∗ = −0.7 and K̄ = 10,000, while k and R are used for normalization according to Sec. 4.6
and remain unspecified. Further, stabilization scheme ‘A-s’ is used with µ̄ = 0.01. The FE
28

Figure 11: Cell budding: initial configuration, FE discretization and boundary conditions. The
boundary normal is fixed and the boundary nodes are only free to move in the radial direction.

discretization shown in Fig. 11, consisting of five NURBS patches, is used. Where the patches
meet, constraint (79) is added to ensure rotational continuity and moment transfer. Constraint
(79) is also used to fix the boundary normal. The actual FE mesh is much finer than in Fig. 11
and uses 12228 elements (64 times more than in the figure). The penalty parameter of the
rotational constraint is ¯ = 6,400. Gaussian quadrature is considered, using 3 × 3 points for
surface integrals and 4 points for line integrals.
5.4.3

Bud shapes

In past numerical studies, axisymmetric bud shapes have been reported, e.g. Walani et al. (2015).
These shapes should be a natural solution due to the axisymmetry of the problem. However, as
is shown below, non-axisymmetric solutions are also possible, and in fact energetically favorable,
indicating that axisymmetric solutions can become unfavored. This is illustrated by considering
case
1
2
3
4
5

bud shape
axisymmetric
general
general
general
general

H0 region
circle
circle
ellipse
ellipse
ellipse

stabilization
A-s
A-s
A-s
A-st
a-st

µ̄
0.01
0.01
0.01
10
1250

in-plane stress state
hydro-static
hydro-static
hydro-static
elastic shear
viscous shear

Table 3: Cell budding: different physical test cases considered.
the five different physical test cases shown in Tab. 3 and discussed in the following:
Case 1: The deformation is constrained to remain axisymmetric (the FE nodes are only allowed
to move in radial direction5 ). The resulting deformation for H̄0 = −25 is shown in Fig. 12 and
in the supplemental movie file bud1.mpg.
Case 2: The deformation is not constrained. Consequently, the non-axisymmetric bud shape
shown in Fig. 13 appears. Viewed from the top, the bud takes the shape of a ‘+’. Essentially,
the initially circular region, where H0 is prescribed, tries to flow away in order to lower the
elastic energy. The flow leads to large mesh distortions and consequently the simulation crashes
at H̄0 = −23.17. It is reasonable to suspect that the ‘+’ shape is a consequence of the particular
discretization. To confirm this, the following case is considered:
5

For the given mesh, this does not enforce axisymmetry exactly as a close inspection of the results shows.

29

Figure 12: Cell budding: Axisymmetric case at H̄0 = −5, −10, −15, −20, −25 (left to right):
3D and side view of deformation and curvature H̄. Here H̄ ∈ [−15.0, 0.31].

Figure 13: Cell budding: Unconstrained, perfect case at H̄0 = −5, −10, −15, −20, −23.17 (left
to right): 3D and side view of deformation and curvature H̄. Here H̄ ∈ [−24.5, 6.12].

Case 3: Instead of a circle, H0 is prescribed within an imperfect circle, i.e. an ellipse with
half-axes a = 0.22R and b = 0.18R. The bud now flows into a ‘−’ shape, as is shown in Fig. 14.
The distortion of the mesh seen on the right clearly shows how the material flows outward in
horizontal direction and back inward in vertical direction. The mesh distortion becomes so
large that the simulation crashes at H̄0 = −11.98. In the current case the flow is not resisted
mechanically (since µ is so low). Resistance is provided either by shear stiffness (e.g. due to an
underlying cytoskeleton) or by viscosity. This is considered in the remaining two cases.
Case 4: Fig. 15 shows the deformation for the same case as before (a = 0.22R and b = 0.18R)
considering now a shear resisting membrane. Model ‘A-st’6 is now considered with µ̄ = 10. The
shear resistance prevents the unbounded spreading of the bud observed earlier. Now the bud
remains localized in the center. The bud starts growing in an almost circular fashion, but then
degenerates into the shape shown in Fig. 15. The process suggests that the initial imperfection is
only a trigger for the final shape, but does not affect it in magnitude. The evolving deformation
for case 4 is also shown in the supplemental movie file bud4.mpg.
Case 5: The final case also considers the imperfect circle from before (a = 0.22R and b =
6

The shear stresses are now physical and need to be applied both in-plane and out-of-plane.

30

Figure 14: Cell budding: Unconstrained, imperfect case at H̄0 = −4, −8, −11.98 (left to right):
3D and side view of deformation and curvature H̄. Here H̄ ∈ [−12.5, −0.27].

Figure 15: Cell budding: Shear resistant, imperfect case at H̄0 = −5, −10, −15, −20, −25 (left
to right): 3D and side view of deformation and curvature H̄. Here H̄ ∈ [−22.8, 4.61].

0.18R), but now provides significant shear stress through physical viscosity. This is captured
through model ‘a-st’ using relation (56) with µ̄ = 1250 and a load stepping increment for H0
of ∆H̄0 = 0.02 (such that ν = 25k/L3 /Ḣ0 , where Ḣ0 is the rate with which the spontaneous
curvature is prescribed). The evolution of the bud with H0 is shown in Fig. 16 and in the
supplemental movie file bud5.mpg. Again, the bud starts growing in an almost circular fashion,
but then degenerates into a plate-like shape. If H0 is kept fixed over time, the solution of case
5 will relax to the solution of case 3.
5.4.4

Surface energy

By examining the surface energy
Z
Π :=

W dA ,

(104)

S0

it can be seen that the non-axisymmetric shapes are preferred. As Fig. 17 shows case 2 and
3 have much lower surface energy than case 1. The difference becomes especially large below
H0 = 4/R, when the deformations become large. As the system tries to minimize energy, this
31

Figure 16: Cell budding: Viscous, imperfect case at H̄0 = −5, −10, −15, −20, −25 (left to right):
3D and side view of deformation and curvature H̄. Here H̄ ∈ [−24.4, 4.29].
9

9
axisymmetric
unconstrained
imperfect

8

7
6
5
4
3

6
5
4
3

2

2

1

1

0
−15

a.

axisymmetric
unconstrained
imperfect

7
elastic energy Π [ k ]

bending energy Πb [ k ]

8

−10
−5
prescribed spontaneous curvature H0 [ 1/R ]

0
−15

0

b.

−10
−5
prescribed spontaneous curvature H0 [ 1/R ]

0

Figure 17: Cell budding: elastic surface energy vs. spontaneous curvature for the three cases
shown in Figs. 12, 13 and 14: a. bending energy; b. total energy (bending + areal part).

shows that the axisymmetric bud shape of case 1 is not a preferred solution. As the figure
shows, almost all of the energy goes into bending (contribution Jw in (27)), as the areal part
(contribution Kg 2 /2 in (27)) becomes negligible for near incompressibility.
5.4.5

Surface tension

One of the advantages of the proposed finite element formulation is that the surface tension γ
can be studied. This is done here for cases 1, 4 and 5 listed above. It was shown before that
the surface tension is not uniform under protein-induced spontaneous curvature through FE
simulations utilizing Monge patches (Rangamani et al., 2014). The current simulations confirm
this result and in addition yield further understanding for large deformations.
First, the axisymmetric case (case 1) is examined in Fig. 18. As seen γ becomes maximum
within the protein patch. This maximum is constant across the patch for low H0 . This changes
for increasing H0 , where a distinct maximum appears in the center.7 At a certain level of H0
7

It can be seen that the distribution of γ is not exactly axisymmetric. This is due to the inexact enforcement

32

Figure 18: Cell budding: Axisymmetric case at H̄0 = −5, −10, −15, −20, −25 (left to right):
3D and side view of deformation and surface tension γ̄.

rupture would occur here, depending on the strength of the lipid bilayer.
Second, the shear stiff case (case 4) is examined in Fig. 19. As the plate-like bud appears, the

Figure 19: Cell budding: Shear resistant, imperfect case at H̄0 = −5, −10, −15, −20, −25 (left
to right): 3D and side view of deformation and surface tension γ̄.

maximum of γ moves away from the center and is concentrated at the end of the plate. The
behavior is similar for the viscous case (case 5) shown in Fig. 20. Both case 4 and 5 show
that the extrema of γ are concentrated in very small regions associated with large curvatures.
However, these peak values are still much lower than the more distributed peak values of case 1.
The lower values are not surprising, as the system has much lower energy for non-axisymmetric
shapes. For the non-axisymmetric cases with shear stiffness (case 4) and viscosity (case 5), the
resulting surface tensions are of similar magnitude. Only the shapes are different.
5.4.6

Effective shear stiffness

Fig. 21 shows the sign of the effective shear stiffness µeff for cases 1 and 4.8 It can be seen that for
of axisymmetry noted in footnote 4.
8
For case 4, the physical shear stiffness µ̄ = 10 needs to be added to Eq. (39).

33

Figure 20: Cell budding: Viscous, imperfect case at H̄0 = −5, −10, −15, −20, −25 (left to right):
3D and side view of deformation and surface tension γ̄.

Figure 21: Cell budding: sign of µ̄eff for cases 1 (top) and 4 (bottom) at H̄0 =
−5, −10, −15, −20, −25 (left to right).

both cases, µeff is negative in some regions. The minimum values of µ̄eff reaches down to −150 in
case 1 and to −50 in case 4. Since the simulations run stably, even though µeff < 0, there must
be another stabilizing effect. It is probable that this related to the geometry: The figure shows
that the regions of negative µeff are all convex. The condition µeff < 0 is therefore not sufficient
for unstable behavior. Since stabilization may be provided naturally, the computations can
sometimes be performed with no stabilization scheme. In practice it is however recommended
to add one of the stabilization schemes proposed in Sec. 3. For example, when both H and κ
are zero neither the Helfrich model nor the geometry can be expected to stabilize the structure.
The stabilization parameters (µ and nt ) can be picked such that the stabilization scheme does
not affect the physical behavior. This is the case for the stabilization chosen here.
5.4.7

Influence of the area-compressibility

As a final study, we investigate the influence of the area-compressibility (parameter K) and
compare the behavior of model (27), that depends on K, with (29), where K is infinite. For
this purpose, model (29) is discretized according to Sec. 4.3. In theory, the behavior of model

34

(27), should approach that of model (29) as K → ∞.
For case 1, there is no significant difference between the two models as Fig. 22 shows. Only

Figure 22: Cell budding: influence of K for the axisymmetric case (case 1) for H̄0 = −15
(top) and H̄0 = −25 (bottom) considering K̄ = 103 , 104 , 105 , 106 , ∞ (left to right). The color
shows γ̄.

for the lowest K, differences appear. Between all other cases, the bud shape only changes
minimally. However, increasing K in formulation (27) leads to oscillations in γ.
For case 5, a strong dependency on K appears, as Fig. 23 shows. Now, the deformation clearly

Figure 23: Cell budding: influence of K for the viscous case (case 5) for H̄0 = −15 (top) and
H̄0 = −25 (bottom) considering K̄ = 103 , 104 , 105 , 106 , ∞ (left to right). The color shows γ̄.
does not converge with increasing K in (27). Instead, the solution for K̄ = 103 is closest to
the solution from model (29). The figure also shows that oscillations appear in the solution
from model (27) as K increases. Essentially, the problem of model (27) is that even though
J converges to 1, as Fig. 24 shows, the pressure q does not converge. This can be seen in
Fig. 25 for the pure bending problem of Sec. 5.1, which has an analytical solution for q. As the
figure shows, there is an optimal value of K, where the pressure error is minimal and beyond
which it diverges. This is different to examples in computational contact mechanics, where
both kinematic and kinetic variables can converge with increasing penalty parameter(Sauer and
De Lorenzis, 2015). So even though model (27) is simpler and more efficient, it has to be used
35

Figure 24: Cell budding: influence of K for case 5 for H̄0 = −25 considering K̄ =
103 , 104 , 106 , ∞ (left to right). The color shows the relative error in J given by |J − 1|.

4

6

10

10
PM: m=1×4
PM: m=1×8
PM: m=1×32
PM: m=1×64
LM: m=1×4
LM: m=1×8
LM: m=1×32
LM: m=1×64

max |1−J|

0

10

4

Max rel. error of pressure

2

10

−2

10

−4

2

10

0

10

−2

10

−4

10

10
0

10

a.

10

2

4

6

0

10
10
10
penalty parameter K [2k/L2]

10

b.

2

4

6

10
10
10
penalty parameter K [2k/L2]

Figure 25: Pure bending: accuracy of the penalty regularization (27) compared to the Lagrange
multiplier formulation (29): a. error in stretch J; b. error in pressure q. As seen, J can be more
accurate, while q does not converge with K. There is rather an optimum K that minimizes the
pressure error.

with care, as the solution can become wrong. If an implementation of (29) is not available, the
mesh convergence of model (27) should be examined

6

Conclusion

This paper presents a general computational formulation for lipid bilayers based on general thinshell kinematics, the Helfrich bending model, and 3D, LBB-conforming, C 1 -continuous NURBSbased finite elements. The rotational continuity of the formulation is ensured by using a rotational constraint across the patch boundaries in the FE mesh. Two cases are considered in order
to model the in-plane membrane response: area-compressibility and area-incompressibility, and
for both suitable FE formulations are presented. Since the formulation lacks shear stiffness,
several shear stabilization schemes are proposed for quasi-static computations. They are based
on adding either numerical stiffness (class A), numerical viscosity (class a) or performing a
projection of the solution onto the shell manifold (class P). The numerical viscosity scheme
can also be used to model physical viscosity as the example in Fig. 16 illustrates, while the
numerical stiffness scheme can also be used to model physical stiffness as the example in Fig. 15
illustrates. It is further shown that the Helfrich bending model provides intrinsic shear stiffness
as long as the surface curvature is non-zero. Altogether, four different computational examples
are considered in order to verify the formulation and to study its physical behavior. The last ex36

ample shows that the 3D budding behavior of lipid bilayers – as described by the Helfrich model
– can become very complex, even though the model is purely mechanical and does not account
for other effects, such as diffusion or temperature. Computational challenges arise especially for
area-incompressible, non-axisymmetric bilayer shapes. For these a penalty regularization can
give misleading results.
The paper shows that the proposed computational model is not only capable of capturing
complex deformations, but is also a suitable tool to analyse and understand the stresses and
energetics of lipid bilayers. Still, more work is needed in order to further advance the computational modeling of lipid bilayers. For example, the modeling of surface flow, protein diffusion,
multicomponent lipids and thermal fluctuations would be useful extensions of the present work.

Acknowledgements
The authors are grateful to the German Research Foundation (DFG) for supporting this research
under grants GSC 111 and SA1822/5-1. They acknowledge support from the University of
California at Berkeley and from Director, Office of Science, Office of Basic Energy Sciences,
Chemical Sciences Division, of the U.S. Department of Energy under contract No. DE AC0205CH11231. Further, they thank the graduate student Yannick Omar for checking the theory.

A

e
Linearization of finto

e
Alternatively to Eq. (69.2), finto
can be also written as
Z
e
finto = −
τ αβ NT (n ⊗ n) N,αβ dA xe .
Ωe0

(105)

e
thus becomes
The linearization of finto
e (x + ∆x ) ≈ f e (x ) + ∆f e
finto
e
e
into e
into ,

with
e
∆finto

Z
= −
Z
−
Ωe0

Z

∆τ αβ bαβ NT n dA
N (n ⊗ n) N,αβ dA ∆xe −
e
Ω0
Z
αβ
T
τ N (n ⊗ ∆n) aα,β dA −
τ αβ bαβ NT ∆n dA .
τ

Ωe0

αβ

(106)

T

(107)

Ωe0

Inserting
∆τ αβ = cαβγδ aγ · N,δ ∆xe

(108)

∆n = −aγ (n · ∆aγ )

(109)

e
∆finto
= keinto ∆xe ,

(110)

(Sauer et al., 2014) and
(Wriggers, 2006), we get
with the tangent matrix
Z
Z
αβ
T
e
kinto = −
τ N (n ⊗ n) N,αβ dA −
cαβγδ bαβ NT (n ⊗ aγ ) N,δ dA
Ωe0
Ωe0
Z
Z
+
τ αβ Γγαβ NT (n ⊗ n) N,γ dA +
τ αβ bαβ NT (aγ ⊗ n) N,γ dA .
Ωe0

Ωe0

37

(111)

Here, Γγαβ := aγ · aα,β defines the Christoffel symbol of the second kind. The tensor components cαβγδ are given in Sauer and Duong (2015) for various material models. For an efficient
implementation, the contractions appearing above should be worked out analytically.

B

Linearization of ge

The vector ge is independent of the Lagrange multiplier and thus only depends on xe . The
linearization of ge thus is
ge (xe + ∆xe ) ≈ ge (xe ) + ∆ge ,
(112)
with

Z

∆ge

=

LT ∆g dA ,

(113)

J αβ
a ∆aαβ .
2

(114)

Ωe0

and (Sauer and Duong, 2015)
∆g = ∆J =

Inserting the discretisation of ∆aαβ (Duong et al., 2015), and exploiting the symmetry of aαβ ,
we get the approximation
∆g ≈ J aα · N,α ∆xe ,
(115)
such that
∆ge = keg ∆xe ,

(116)

with the tangent matrix
keg

Z
:=
Ωe0

LT aα · N,α J dA .

(117)

References
Agrawal, A. and Steigmann, D. (2009). Modeling protein-mediated morphology in biomembranes. Biomech. Model. Mechanobiol., 8(5):371–379.
Aris, R. (1989). Vectors, tensors and the basic equations of fluid mechanics. Dover, Mineola.
Arroyo, M. and DeSimone, A. (2009). Relaxation dynamics of fluid membranes. Phys. Rev. E,
79:031915.
Babuška, I. (1973). The finite element method with Lagrangian multipliers. Num. Math.,
20:179–192.
Baesu, E., Rudd, R. E., Belak, J., and McElfresh, M. (2004). Continuum modeling of cell
membranes. Int. J. Non-lin. Mech., 39:369–377.
Bathe, K.-J. (1996). Finite Element Procedures. Prentice-Hall, New Jersey.
Borden, M. J., Scott, M. A., Evans, J. A., and Hughes, T. J. R. (2011). Isogeometric finite
element data structures based on Bezier extraction of NURBS. Int. J. Numer. Meth. Engng.,
87:15–47.
Budin, I., Bruckner, R. J., and Szostak, J. W. (2009). Formation of protocell-like vesicles in a
thermal diffusion column. Journal of the American Chemical Society, 131(28):9628–9629.
38

Buser, C. and Drubin, D. G. (2013). Ultrastructural imaging of endocytic sites in saccharomyces cerevisiae by transmission electron microscopy and immunolabeling. Microscopy and
Microanalysis, 19(02):381–392.
Canham, P. B. (1970). The minimum energy of bending as a possible explanation of the
biconcave shape of the human red blood cell. J. Theoret. Biol., 26:61–81.
Cirak, F. and Ortiz, M. (2001). Fully C1 -conforming subdivision elements for finite elementdeformation thin-shell analysis. Int. J. Numer. Meth. Engng, 51:813–833.
Cottrell, J. A., Hughes, T. J. R., and Bazilevs, Y. (2009). Isogeometric Analysis. Wiley.
Dao, M., Lim, C. T., and Suresh, S. (2003). Mechanics of the human red blood cell deformed
by optical tweezers. J. Mech. Phys. Solids, 51:2259–2280.
Derényi, I., Jülicher, F., and Prost, J. (2002). Formation and interaction of membrane tubes.
Phy. Rev. Lett., 88(23):238101.
Du, Q. and Wang, X. Q. (2007). Convergence of numerical approximations to a phase field
bending elasticity model of membrane deformations. Int. J. Numer. Anal. Model., 4(3-4):441–
459.
Duong, T. X., Roohbakhshan, F., and Sauer, R. A. (2015). A new rotation-free isogeometric thin
shell formulation and a corresponding continuity constraint for patch boundaries. submitted.
Elliott, C. M. and Stinner, B. (2010). Modeling and computation of two phase geometric
biomembranes using surface finite elements. J. Comp. Phys., 229(18):6585–6612.
Fawcett, D. W. (1981). The Cell. Saunders, Philadelphia.
Feng, F. and Klug, W. S. (2006). Finite element modeling of lipid bilayer membranes. J.
Comput. Phys., 220:394–408.
Griparic, L. and van der Bliek, A. M. (2001). The many shapes of mitochondrial membranes.
Traffic, 2(4):235–244.
Gruenberg, J. and Stenmark, H. (2004). The biogenesis of multivesicular endosomes. Nature
Reviews Molecular Cell Biology, 5(4):317–323.
Gu, R., Wang, X., and Gunzburger, M. (2014). Simulating vesicle-substrate adhesion using two
phase field functions. J. Comp. Phys., 275:626–641.
Helfrich, W. (1973). Elastic properties of lipid bilayers: Theory and possible experiments. Z.
Naturforsch., 28c:693–703.
Hu, J., Shibata, Y., Zhu, P.-P., Voss, C., Rismanchi, N., Prinz, W. A., Rapoport, T. A., and
Blackstone, C. (2009). A class of dynamin-like gtpases involved in the generation of the
tubular er network. Cell, 138(3):549–561.
Itoh, T., Erdmann, K. S., Roux, A., Habermann, B., Werner, H., and De Camilli, P. (2005).
Dynamin and the actin cytoskeleton cooperatively regulate plasma membrane invagination
by bar and f-bar proteins. Developmental cell, 9(6):791–804.
Jarić, M., Seifert, U., Wirtz, W., and Wortis, M. (1995). Vesicular instabilities: The prolateto-oblate transition and other shape instabilities of fluid bilayer membranes. Phys. Rev. E,
52(6):6623–6634.

39

Jie, Y., Quanhui, L., Jixing, L., and Zhong-Can, O.-Y. (1998). Numerical observation of nonaxisymmetric vesicles in fluid membranes. Phys. Rev. E, 58(4):4730–4736.
Kahraman, O., Stoop, N., and Müller, M. M. (2012). Fluid membrane vesicles in confinement.
New J. Phys., 14:095021.
Karotki, L., Huiskonen, J. T., Stefanand, C. J., Ziólkowska, N. E., Roth, R., Surma, M. A.,
Krogan, N. J., Emr, S. D., Heuser, J., Grünewald, K., and Walther, T. C. (2011). Eisosome
proteins assemble into a membrane scaffold. J. Cell Biol., 195(5):889–902.
Keren, K. (2011). Cell motility: the integrating role of the plasma membrane. European
Biophysics Journal, 40(9):1013–1027.
Kishimoto, T., Sun, Y., Buser, C., Liu, J., Michelot, A., and Drubin, D. G. (2011). Determinants
of endocytic membrane geometry, stability, and scission. Proc. Natl. Acad. Sci., 108:E979–
988.
Kloeppel, T. and Wall, W. A. (2011). A novel two-layer, coupled finite element approach for
modeling the nonlinear elastic and viscoelastic behavior of human erythrocytes. Biomech.
Model. Mechanobiol., 10(4):445–459.
Koster, G., VanDuijn, M., Hofs, B., and Dogterom, M. (2003). Membrane tube formation from
giant vesicles by dynamic association of motor proteins. Proceedings of the National Academy
of Sciences, 100(26):15583–15588.
Kozlov, M. M., Campelo, F., Liska, N., Chernomordik, L. V., Marrink, S. J., and McMahon,
H. T. (2014). Mechanisms shaping cell membranes. Current opinion in cell biology, 29:53–60.
Kukulski, W., Schorb, M., Kaksonen, M., and Briggs, J. A. G. (2012). Plasma membrane reshaping during endocytosis is revealed by time-resolved electron tomography. Cell, 150(3):508–
520.
Lau, C., Brownell, W. E., and Spector, A. A. (2012). Internal forces, tension and energy density
in tethered cellular membranes. J. Biomech., 45(7):1328–1331.
Li, H. and Lykotrafitis, G. (2012). Two-component coarse-grained molecular-dynamics model
for the human erythrocyte membrane. Biophys. J., 102(1):75–84.
Ma, L. and Klug, W. S. (2008). Viscous regularization and r-adaptive meshing for finite element
analysis of lipid membrane mechanics. J. Comput. Phys., 227:5816–5835.
McMahon, H. T. and Gallop, J. L. (2005). Membrane curvature and mechanisms of dynamic
cell membrane remodelling. Nature, 438(7068):590–596.
Peng, Z., Asaro, R. J., and Zhu, Q. (2010). Multiscale simulation of erythrocyte membranes.
Phys. Rev. E, 81:031904.
Peter, B. J., Kent, H. M., Mills, I. G., Vallis, Y., Butler, P. J. G., Evans, P. R., and McMahon,
H. T. (2004). Bar domains as sensors of membrane curvature: the amphiphysin bar structure.
Science, 303(5657):495–499.
Rahimi, M. and Arroyo, M. (2012). Shape dynamics, lipid hydrodynamics, and the complex
viscoelasticity of bilayer membranes. Phys. Rev. E, 86:011932.
Ramakrishnan, N., Kumar, P. B. S., and Ipsen, J. H. (2010). Monte carlo simulations of fluid
vesicles with in-plane orientational ordering. Phys. Rev. E, 81:041922.
40

Rangamani, P., Agrawal, A., Mandadapu, K. K., Oster, G., and Steigmann, D. J. (2013). Interaction between surface shape and intra-surface viscous flow on lipid membranes. Biomech.
Model. Mechanobiol., 12(4):833–845.
Rangamani, P., Mandadapu, K. K., and Oster, G. (2014). Protein-induced membrane curvature
alters local membrane tension. Biophysical journal, 107(3):751–762.
Rangarajan, R. and Gao, H. (2015). A finite element method to compute three-dimensional
equilibrium configurations of fluid membranes: Optimal parameterization, variational formulation and applications. J. Comp. Phys., 297:266–294.
Rim, J. E., Purohit, P. K., and Klug, W. S. (2014). Mechanical collapse of confined fluid
membrane vesicles. Biomech. Model. Mechanobio., 13(6):1277–1288.
Rosolen, A., Peco, C., and Arroyo, M. (2013). An adaptive meshfree method for phase-field
models of biomembranes. Part I: Approximation with maximum-entropy basis functions. J.
Comput. Phys., 249:303–319.
Salac, D. and Miksis, M. (2011). A level set projection model of lipid vesicles in general flows.
J. Comput. Phys., 230(22):8192–8215.
Sauer, R. A. (2014). Stabilized finite element formulations for liquid membranes and their
application to droplet contact. Int. J. Numer. Meth. Fluids, 75(7):519–545.
Sauer, R. A. and De Lorenzis, L. (2015). An unbiased computational contact formulation for
3D friction. Int. J. Numer. Meth. Engrg., 101(4):251–280.
Sauer, R. A. and Duong, T. X. (2015). On the theoretical foundations of solid and liquid shells.
Math. Mech. Solids, published online, DOI: 10.1177/1081286515594656.
Sauer, R. A., Duong, T. X., and Corbett, C. J. (2014). A computational formulation for solid and
liquid membranes based on curvilinear coordinates and isogeometric finite elements. Comput.
Methods Appl. Mech. Engrg., 271:48–68.
Shemesh, T., Klemm, R. W., Romano, F. B., Wang, S., Vaughan, J., Zhuang, X., Tukachinsky,
H., Kozlov, M. M., and Rapoport, T. A. (2014). A model for the generation and interconversion of er morphologies. Proceedings of the National Academy of Sciences, 111(49):E5243–
E5251.
Shi, Z. and Baumgart, T. (2015). Membrane tension and peripheral protein density mediate
membrane shape transitions. Nature communications, 6:5974.
Shibata, Y., Hu, J., Kozlov, M. M., and Rapoport, T. A. (2009). Mechanisms shaping the
membranes of cellular organelles. Annual Review of Cell and Developmental Biology, 25:329–
354.
Shibata, Y., Voeltz, G. K., and Rapoport, T. A. (2006). Rough sheets and smooth tubules.
Cell, 126(3):435–439.
Steigmann, D. J. (1999). Fluid films with curvature elasticity. Arch. Rat. Mech. Anal., 150:127–
152.
Tasso, I. V. and Buscaglia, G. C. (2013). A finite element method for viscous membranes.
Comput. Meth. Appl. Mech. Engrg., 255:226–237.
Terasaki, M., Chen, L. B., and Fujiwara, K. (1986). Microtubules and the endoplasmic reticulum
are highly interdependent structures. The Journal of Cell Biology, 103(4):1557–1568.
41

Walani, N., Torres, J., and Agrawal, A. (2015). Endocytic proteins drive vesicle growth via
instability in high membrane tension environment. Proceedings of the National Academy of
Sciences, 112(12):E1423–E1432.
Walther, T. C., Brickner, J. H., Aguilara, P. S., Bernales, S., Pantoja, C., and Walter, P. (2006).
Eisosomes mark static sites of endocytosis. Nature, 439:998–1003.
Waterman-Storer, C. M. and Salmon, E. D. (1998). Endoplasmic reticulum membrane tubules
are distributed by microtubules in living cells using three distinct mechanisms. Current
Biology, 8(14):798–807.
Wriggers, P. (2006). Computational Contact Mechanics. Springer, 2nd edition.
Wriggers, P. (2008). Nonlinear Finite Element Methods. Springer.
Zimmerberg, J. and Kozlov, M. M. (2006). How proteins produce cellular membrane curvature.
Nature Reviews Molecular Cell Biology, 7(1):9–19.

42

