# Learning Stable Koopman Embeddings

Fletcher Fan, Bowen Yi, David Rye, Guodong Shi and Ian R. Manchester This work
was supported by the Australian Research Council.The authors are with the
Australian Center for Field Robotics and Sydney Institute for Robotics and
Intelligent Systems, The University of Sydney, NSW 2006, Australia.
Corresponding author’s email<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a new data-driven method for learning stable models
of nonlinear systems. Our model lifts the original state space to a higher-
dimensional linear manifold using Koopman embeddings. Interestingly, we prove
that every discrete-time nonlinear contracting model can be learnt in our
framework. Another significant merit of the proposed approach is that it
allows for _unconstrained_ optimization over the Koopman embedding and
operator jointly while enforcing stability of the model, via a direct
parameterization of stable linear systems, greatly simplifying the
computations involved. We validate our method on a simulated system and
analyze the advantages of our parameterization compared to alternatives.

##  I Introduction

The problem of fitting models to data generated from dynamical systems, known
as system identification, is ubiquitous in science and engineering. One
important consideration in system identification is the stability of the
model. In many applications, the fitted model is used for prediction of future
behaviors of the system, and an unstable model would erroneously produce
unbounded predictions.

There are many different forms of stability for nonlinear systems. In this
paper, we consider contraction, also known as incremental stability, which can
be viewed as a “strong” type of stability for dynamical systems that studies
the convergence between _any_ two trajectories of the given system [10]. There
has been much prior work on learning contracting models in system
identification for various model classes, including polynomial models [26,
27], Gaussian mixture models [20] and neural network models [21, 22, 13]. In
this paper, we propose a new class of contracting nonlinear models that
combines the expressiveness of neural networks with the strong stability
guarantees associated with linear systems. Furthermore, we propose a learning
framework that fits this class of models to data via an unconstrained
optimization problem.

Our work bridges the gap between learning stable nonlinear models and
approximating the Koopman operator [7], an infinite-dimensional linear
operator that can describe the dynamics of any nonlinear system by embedding
it in a higher-dimensional space. There has been growing interest in data-
driven methods that estimate finite-dimensional approximations of the Koopman
operator and its eigenfunctions [24, 28, 4], motivated by the appeal of being
able to apply linear systems analysis to complex nonlinear systems. Due to
this merit, Koopman-based methods have been developed for system
identification [15], state observation and control [8] of nonlinear systems.

In Koopman identification approaches, a central problem is how to learn the
Koopman embedding from data. Recently many methods [9, 11, 14, 18, 19, 25, 29]
have been proposed to address this problem, however most of them do not
consider the stability of the learned model. Unstable learned models may have
serious robustness issues, particularly when applied to dissipative physical
systems, making them unsuitable for practical use. We aim to address this
issue by imposing stability constraints on the Koopman model.

Our model class is motivated by recent work [30] showing that, for continuous-
time (CT) nonlinear systems, there is an equivalence between the Koopman and
contraction approaches for stability analysis under some mild technical
assumptions. We extend this equivalence result to discrete-time (DT) systems
in this paper, and provide an algorithmic framework for learning Koopman
models with contracting properties.

The main contributions of this paper are threefold:

  * C1

We propose a novel Koopman learning framework that jointly models the Koopman
operator and embedding from data, while imposing the model
stability/contraction constraint.

  * C2

We prove that every nonlinear discrete-time contracting model can be learnt in
our framework in an arbitrarily large compact set, which may be viewed as the
extension of [30] from CT to DT; see Theorem 1.

  * C3

Our work builds on the contracting model class identified in [22], which
allows for _unconstrained optimization_ of objective functions, unlike some
existing parameterizations of Koopman operators, e.g. [12]. As a result, it
significantly simplifies the implementation of optimization algorithms for
learning the model parameters.

The rest of the paper is organized as follows. Section II defines the system
identification problem and provides the background on Koopman operator theory
and data-driven methods for estimating the Koopman operator, and also restates
the main result in [30]. Section III extends [30, Theorem 1] to DT systems.
Section IV defines the model set used in our learning framework, and Section V
defines the optimization problem. Section VII provides some numerical
validations of our framework on a handwriting dataset.

Notation. All mappings and functions are assumed sufficiently smooth. Given
$f:\mathbb{R}^{n}\to\mathbb{R}^{m}$, we denote the gradient operator $\nabla
f:=(\partial f/\partial x)^{\top}$. Let ${\cal F}$ be the space of smooth
real-valued scalar functions $\mathbb{R}^{n}\to\mathbb{R}$. We use
$(\cdot)^{\dagger}$ to denote the Moore-Penrose pseudoinverse of a matrix.
$\lambda_{\tt min}(\cdot)$ and $\lambda_{\tt max}(\cdot)$ respectively
represent the smallest and largest eigenvalue of a square matrix. We use
$\lvert\cdot\rvert$ to denote vector norms, i.e. $\lvert\cdot\rvert_{2}$ is
the vector 2-norm. Sometimes, we may simply write $x(t)$ as $x_{t}$.

##  II Problem Definition and Background

In this paper, we consider the identification of a Koopman embedding and
operator for a discrete-time (DT) autonomous state-space system:

$x(t+1)=f(x(t)),$ (1)

where $x\in\mathbb{R}^{n}$ and $t$ is the timestep. We assume the system (1)
has a single equilibrium at $x_{\star}$, i.e. $f(x_{\star})=x_{\star}$.
Further, we assume the dynamics $f(x)$ are unknown, but we have access to
full-state trajectory data $\\{\tilde{x}_{t}\\}^{T}_{t=0}$, generated by
system (1). We are concerned with learning a function $\phi(x)$ (i.e. the
Koopman embedding) that smoothly maps from the original state space
$\mathbb{R}^{n}$ to a possibly higher-dimensional space $\mathbb{R}^{N}$
($N\geq n$), as well as a linear matrix ${\cal A}\in\mathbb{R}^{N\times N}$
(i.e. a finite-dimensional approximation of the Koopman operator) that
describes the evolution of $\phi(x)$ over time.

The Koopman embedding and the matrix ${\cal A}$ in fact form a predictive
model of the system (1), which we define as a Koopman model.

######  Definition 1 (DT Koopman model)

Given a Koopman embedding $\phi(x)$ and matrix ${\cal A}$, the corresponding
Koopman model is:

$x(t)=a(x_{0},t)=\phi^{L}({\cal A}^{t}\phi(x_{0})),$ (2)

where $\phi^{L}:\mathbb{R}^{N}\to\mathbb{R}^{n}$ is a left-inverse of
$\phi(x)$ such that $\phi^{L}(\phi(x))=x$, and $x_{0}$ is an initial
condition.

The problem of learning $\phi(x)$ and ${\cal A}$ can be treated as a
minimization of the prediction error of the Koopman model on the given data
$\\{\tilde{x}_{t}\\}^{T}_{t=0}$.

###  II-A Koopman Operator Theory

Before presenting our theoretical contributions, we provide some background on
the Koopman operator. The Koopman operator was proposed in [7] for CT
dynamical models. First, let us recall the definition of its variant for DT
systems.

###### Definition 2

(Koopman operator) For the DT dynamical model (1), the Koopman operator ${\cal
K}:{\cal F}\to{\cal F}$ is defined by

${\cal K}[\varphi(x)]:=\varphi\circ f(x)$ (3)

for $\varphi\in{\cal F}$, assuming that the system has a unique solution
$\forall t\in\mathbb{N}$. We term the scalar real-valued function
$\varphi:\mathbb{R}^{n}\to\mathbb{R}$ an _observable_.

Since the Koopman operator is defined on the functional space, it is infinite-
dimensional. It is also easy to verify that the Koopman operator is linear,
i.e. ${\cal K}[k_{1}\varphi_{1}+k_{2}\varphi_{2}]=k_{1}{\cal
K}[\varphi_{1}]+k_{2}{\cal K}[\varphi_{2}]$ for any $k_{1},k_{2}\in\mathbb{R}$
and $\varphi_{1},\varphi_{2}\in{\cal F}$. This property makes Koopman methods
widely popular in the analysis of dynamical models. Despite the infinite
dimension of the Koopman operator, some key properties of a given nonlinear
dynamical model—e.g. stability and dynamical behaviours—can be captured by a
few particular functions, i.e. the Koopman eigenfunctions.

###### Definition 3

(Koopman eigenfunction) A Koopman eigenfunction is a non-zero observable
$\phi_{\lambda}\in{\cal F}/\\{0\\}$ satisfying

${\cal K}[\phi_{\lambda}(x)]=\lambda\phi_{\lambda}(x)$ (4)

for some $\lambda\in\mathbb{C}$, which is the associated Koopman eigenvalue.

A Koopman eigenfunction defines a coordinate in which the system trajectories
behave as a linear system. To be precise, define a coordinate
$z_{\lambda}=\phi_{\lambda}(x)$, the dynamics of which are given by

$\displaystyle z_{\lambda}(t+1)=\lambda z_{\lambda}(t),$

with the initial condition $z_{\lambda}(0)=\phi_{\lambda}(x(0))$. Indeed, the
definition (4) is equivalent to solving the algebraic equation

$\phi_{\lambda}(f(x))=\lambda\phi_{\lambda}(x),\quad\forall
x\in\mathbb{R}^{n}$

if the DT dynamical model (1) is prior.

###  II-B Dynamic Mode Decomposition

In the problem of system identification, we are more interested in finding the
Koopman eigenfunctions and eigenvalues only from the collected data set
$\\{\tilde{x}_{t}\\}_{t=0}^{T}$, for which dynamic mode decomposition (DMD)
provides an efficient data-driven approach to approximating the Koopman
operator [24].

In DMD, usually some heuristically _predetermined_ , sufficiently rich
observables $\phi_{1},\ldots\phi_{N}$ ($N\gg n$)—rather than Koopman
eigenfunctions—are involved to learn the nonlinearity in the dynamical model.
The task in the DMD method is to seek a matrix ${\cal A}\in\mathbb{R}^{N\times
N}$ in order to obtain a finite-dimensional approximation of ${\cal K}$, which
minimizes the following:

$\sum_{j=0}^{T}|\phi(x(t+1))-{\cal A}\phi(x(t)))|^{2}_{2},$ (5)

in which we have defined $\psi:=[\phi_{1},\ldots,\phi_{N}]^{\top}$. The least
square problem (5) has a unique solution

${\cal A}=Y_{1}Y_{2}^{\dagger}$ (6)

with

$\displaystyle Y_{1}$
$\displaystyle:=\begin{bmatrix}\phi(x(1))&\ldots&\phi(x(T))\end{bmatrix}$
$\displaystyle Y_{2}$
$\displaystyle:=\begin{bmatrix}\phi(x(0))&\ldots&\phi(x(T-1))\end{bmatrix}$

if $Y_{2}$ is full row rank. DMD is a simple, efficient method to approximate
the Koopman operator, but two issues arise:

  * 1)

In the DMD method, the observables $\psi$ are predetermined, which
significantly affects the learning accuracy, but in the literature the
selection of observables usually done in a heuristic manner. Since these
observables are closely connected to the Koopman eigenfunctions for a given
dynamical model, a natural question is: can the observables and the matrix
${\cal A}$ be learnt concurrently to improve accuracy?

  * 2)

For a stable dynamical model, the above least square solution may yield an
unstable model due to various kinds of perturbations in the data set
$\\{\tilde{x}_{t}\\}^{T}_{t=0}$, which would be unacceptable in many
applications. Hence, imposing stability constraints is an important
consideration in learning algorithms.

The main motivation of the paper is to address the above issues and present a
novel Koopman learning framework.

###  II-C Contraction analysis

In this paper, we are interested in stable nonlinear models. Indeed, there are
many different forms of stability for nonlinear systems; we focus on
contracting systems [10].

Contraction analysis provides another way to study nonlinear systems by means
of linear systems theory _exactly_ and _globally_. In contraction analysis we
are concerned with the differential dynamics of a given system, which is
indeed a linear time-varying (LTV) system. The differential dynamics of the
model (1) are given by

$\delta x(t+1)={\partial f\over\partial x}(x(t))\delta x(t),$ (7)

with $\delta x\in\mathbb{R}^{n}$ representing the infinitesimal displacement.
Informally, if the LTV system (7) is exponentially stable along any feasible
trajectories $x(t)$, we can say the system (1) is contracting. Its formal
definition is given as follows.

###### Definition 4

Given the DT system (1), if there exists a uniformly bounded metric $M(x)$,
i.e. $a_{1}I_{n}\preceq M(x)\preceq a_{2}I_{n}$ for some $a_{2}\geq a_{1}>0$,
guaranteeing

${\partial f\over\partial x}(x(t))^{\top}M(x(t+1)){\partial f\over\partial
x}(x(t))-M(x(t))\preceq-\beta M(x(t)),$ (8)

with $0<\beta<1$, then we say that the given system is contracting.

A central result of contraction analysis is that, for contracting systems, all
trajectories converge exponentially to a single trajectory, i.e., for any two
trajectories $x_{a}$ and $x_{b}$

$|x_{a}(t)-x_{b}(t)|\leq a_{0}\beta^{t}|x_{a}(0)-x_{b}(0)|$

for some $a_{0}>0$.

In this paper, we propose an algorithm to learn dynamical models which are
contracting in the sense of Definition 4. Generally speaking, verifying or
guaranteeing contraction for a nonlinear model is non-trivial. Motivated by
the Koopman approach, we instead consider a transformation of the state space
into a higher-dimensional manifold on which the dynamics are linear. A natural
question that arises is the conservativeness of such an approach for verifying
stability. However, it was recently shown in [30] that the Koopman and
contraction approaches are equivalent to each other for stability analysis
when considering CT dynamical models.

##  III Stability Criterion for Discrete-time Koopman Models

In this section, we extend the main result in [30] to DT systems, i.e. the
Koopman and contraction approaches are equivalent for nonlinear stability
analysis. We will show that the model set proposed here can provide sufficient
degrees of freedom for learning nonlinear DT models.

###### Theorem 1

Consider the system (1). Suppose that there exists a mapping
$\phi:\mathbb{R}^{n}\to\mathbb{R}^{N}$ with $N\geq n$ such that

  * D1

There exists Schur stable matrix ${\cal A}\in\mathbb{R}^{N\times N}$
satisfying

$\phi(f(x))={\cal A}\phi(x).$ (9)

  * D2

$\Phi(x):=\nabla\phi(x)^{\top}$ has full column rank, and
$\Phi(x)^{\top}\Phi(x)$ is uniformly bounded.

Then system (1) is contracting with the contraction metric
$\Phi(x)^{\top}P\Phi(x)$, where $P$ is any positive-definite matrix satisfying
$P-{\cal A}^{\top}P{\cal A}\succ 0$. Conversely, if the system (1) is
contracting with the metric $M(x)\in\mathbb{R}^{n\times n}_{\succ 0}$, and
assuming that $f$ is invertible and its inverse $f^{-1}$ is continuous. Then,
in any invariant compact set ${\cal X}\subset\mathbb{R}^{n}$, there exists a
continuous Koopman mapping $\phi:\mathbb{R}^{n}\to\mathbb{R}^{n}$ verifying D1
and D2.

###### Proof:

($\Rightarrow$) From D2 there exists a matrix $P=P^{\top}\succ 0$ satisfying
the Lyapunov condition

$P-{\cal A}^{\top}P{\cal A}\succ Q,$ (10)

for some constant positive definite matrix $Q\succ 0$ without loss of
generality. We define a new coordinate $z:=\phi(x)$ in which the infinitesimal
displacement $\delta z_{t}\in\mathbb{R}^{N}$ at time $t$ is given by

$\delta z_{t}=\Phi(x_{t})\delta x_{t},$ (11)

where $\delta x_{t}$ is an infinitesimal displacement in the $x$-coordinate.

The DT differential dynamics of $x$ can be written as

$\delta x_{t+1}=F(x_{t})\delta x_{t},$ (12)

where $F(x):=\nabla f(x)^{\top}$.

Similarly, for $z$ we have

$\delta z_{t+1}={\cal A}\Phi(x_{t})\delta x_{t}=\Phi(x_{t+1})F(x_{t})\delta
x_{t},$ (13)

where we have used the relations $z_{t+1}=Az_{t}$ and $z_{t+1}=\phi(f(x_{t}))$
in their differential forms. Hence, we obtain

$\Phi(x_{t+1})F(x_{t})={\cal A}\Phi(x_{t}).$

Due to the full column rank of $\Phi(x)$ and (10), it follows that

$\Phi(x_{t})^{\top}(P-{\cal A}^{\top}P{\cal
A})\Phi(x_{t})\succ\Phi(x_{t})^{\top}Q\Phi(x_{t}).$ (14)

Then, by substituting (13), we have

$\displaystyle~{}\Phi(x_{t})^{\top}P\Phi(x_{t})-F(x_{t})^{\top}\Phi(x_{t+1})^{\top}P\Phi(x_{t+1})F(x_{t})$
(15) $\displaystyle\overset{\eqref{eqn:diff_dyn}}{\succ}$
$\displaystyle~{}\Phi^{\top}Q\Phi\succeq{\lambda_{\tt min}(Q)\over\lambda_{\tt
max}(P)}\Phi^{\top}P\Phi.$

Now since $\Phi$ has full column rank and $P\succ 0$, we have
$M(x):=\Phi^{\top}P\Phi\succ 0$. Substituting into (15):

$M(x_{t})-{F}(x_{t})^{\top}M(x_{t+1}){F}(x_{t})\succ\beta M(x_{t}),$ (16)

with $\beta:={\lambda_{\tt min}(Q)/\lambda_{\tt max}(P)}$. By selecting
$Q=\rho P$ with $\rho\in(0,1)$, we have $\beta\in(0,1)$. This is exactly the
contraction condition for the system (1) with respect to the metric $M$.

($\Leftarrow$) For the given DT system, from directly applying the Banach
fixed-point theorem we conclude that there exists a unique fixed-point
$x_{\star}\in{\cal X}$, i.e. $f(x_{\star})=x_{\star}$.

First, we parameterise the unknown mapping $\phi(x)$ as $\phi(x):=x+T(x),$
with a new mapping $T(x)$ to be searched for. Then, the algebraic equation (9)
becomes $T(f(x))+f(x)={\cal A}x+{\cal A}T(x).$ By fixing ${\cal A}=\nabla
f(x_{\star})^{\top}$, from the contraction assumption, we have
$M(x_{\star})-{\cal A}^{\top}M(x_{\star}){\cal A}\succeq\beta M(x_{\star}),$
thus ${\cal A}$ being Schur stable. It yields

$T(f(x))={\cal A}T(x)+H(x),$ (17)

in which we have defined $H(x):={\cal A}x-f(x).$ We make the key observation
that the algebraic equation (17) exactly coincides with the one in the
formulation of the Kazantzis-Kravaris-Luenberger observer for nonlinear DT
systems [2, Eq. (7)]. In our case, the function $H(x)$ is continuous and,
following [2, Theorem 2], we have a feasible solution to (17) as follows:
111The second assumption in [2] holds true in any backward invariant compact
set. Since contracting systems generally cannot guarantee such invariance, we
may modify the dynamics as $x_{t+1}=\rho(x_{t})f(x_{t})$ with
$\rho(x)=\left\\{\begin{aligned}
1,\quad&\mbox{if~{}~{}}x\in\mathtt{cl}(\mathcal{X})\\\ 0,\quad&\mbox{if
~{}~{}}x\notin\mathcal{X}^{\prime}\end{aligned}\right.$ with ${\cal
X}\subset{\cal X}^{\prime}$, and then continue the analysis.

$T(x)=\sum_{j=0}^{+\infty}{\cal A}^{i}H(X(x,-j+1)),$ (18)

with the definition

$X(x,j)=\underbrace{f\circ f\circ\cdots\circ f}_{j~{}\mbox{\small
times}}(x),\quad X(x,-j)=(f^{-1})^{j}(x)$

for $i\in\mathbb{N}_{+}$.

Although $\phi^{0}(x):=x+T(x)$ with $T$ defined above satisfies D1 in the
entire set ${\cal X}$, the condition D2 may be not true. Hence, we need to
modify the obtained $\phi^{0}(x)$. By considering the evolution of the
trajectories in the $x$\- and $z:=\phi(x)$-coordinates respectively, we have

$z(t_{x})=\phi^{0}(x(t_{x}))=\phi^{0}(X(x,t_{x}))={\cal
A}^{t_{x}}\phi^{0}(x),$

with $t_{x}\in\mathbb{N}_{+}$, thus satisfying $\phi^{0}(x)={\cal
A}^{-t_{x}}\phi^{0}(X(x,t_{x})).$ Then, we modify $\phi^{0}(x)$ into

$\phi(x):={\cal A}^{-t_{x}}[X(x,t_{x})+T(X(x,t_{x}))]$ (19)

with a sufficiently large $t_{x}\in\mathbb{N}_{+}$.

Finally, let us check conditions D1 and D2. For the algebraic condition, we
have

$\displaystyle\phi(f(x))$ $\displaystyle={\cal
A}^{-t_{x}}\phi^{0}(X(f(x),t_{x}))$ $\displaystyle={\cal
A}^{-t_{x}}\phi^{0}(f(X(x,t_{x})))$ $\displaystyle={\cal A}^{-t_{x}}\cdot{\cal
A}\phi^{0}(X(x,t_{x}))$ $\displaystyle={\cal A}\phi(x)$

where we have used the fact

$X(f(x),t_{x})=\underbrace{f\circ f\circ\cdots\circ f}_{(j+1)~{}\mbox{\small
times}}=f(X(x,t_{x}))$

in the second equation. Therefore, $\phi(x)$ defined in (19) satisfies the
algebraic equation (9). Regarding D2, let us study the Jacobian of $\phi(x)$
in (19), which is given by

${\partial\phi\over\partial x}(x)={\cal A}^{-t_{x}}\left[I+{\partial
T\over\partial x}(X(x,t_{x}))\right]{\partial X\over\partial x}(x).$

On the other hand, we have that $\nabla_{x}X$ is full rank and

$H(x_{\star})=0,\quad{\partial H\over\partial x}(x_{\star})=0,$

as a result $\nabla T(x_{\star})=0$. If $t_{x}\in\mathbb{N}_{+}$ is
sufficiently large, the largest singular value of $\nabla T(X(x,t_{x}))$ would
be very small, and then the identity part of $\phi(x)$ will dominate
$\nabla\phi(x)$. Hence, $\phi(x)$ is an injection for a large
$t_{x}\in\mathbb{N}_{+}$. ∎

##  IV Model Set

In this section, we define the model set that we optimize over in our learning
framework. We parameterize both the Koopman observables $\phi(x)$ and the
matrix ${\cal A}$ in our model and train them jointly. To the best of our
knowledge, the joint learning of $\phi$ and ${\cal A}$ with the model
stability constraint—as done in the paper—has not been previously considered
in the literature.

Recall our definition of the Koopman model (2). By Theorem 1, the model (2) is
guaranteed to be contracting if Conditions D1 and D2 are satisfied. In the
following, we propose parameterizations of $\phi(x)$ and ${\cal A}$ that
satisfy these conditions.

###  IV-A Parameterization of observables

We propose to parameterize the observables as:

$\phi(x)=Cx+\varphi(x,\theta_{NN}),$ (20)

where $C=[I_{n},0_{n\times(N-n)}]^{\top}$. The nonlinear part $\varphi(x)$ can
be any differentiable function approximator, parameterized by $\theta_{NN}$.
For brevity, we drop the dependence on $\theta_{NN}$ in our notation. In this
paper, we consider $\varphi(x)$ as a feedforward neural network due to its
scalability, but any differentiable function approximator can be used.

The dimensionality of the observables $N$ is a hyperparameter chosen by the
user. For $N=n$, the observables will be of the same form as the constructive
mapping $\phi^{0}(x)=x+T(x)$ in Theorem 1.

In order to reconstruct the original state $x$ from the observables, we need
to train a separate function $\phi^{L}(z)$ to compute the left-inverse of
$\phi(x)$. Indeed, the left invertibility of $\phi$ is necessary for condition
D2. We propose to simply parameterize this left inverse function as another
neural network $\phi^{L}(z,\theta_{L})$.

###### Remark 1

There are many possible parameterizations of the observables that are
compatible with our framework, with Equation (20) being just the one chosen to
mimic the constructive mapping from Theorem 1. For some parameterizations, the
left inverse may be computed analytically and does not have to be modelled as
a separate function. For example, if
$\phi(x)=[x^{\top},\varphi(x)^{\top}]^{\top}$, then the left inverse is simply
$x=C\phi(x)$, where $C=[I,0]$.

###  IV-B Parameterization of the Koopman operator

The Koopman matrix ${\cal A}$ has to satisfy Condition D1 of Theorem 1, i.e.
it must be Schur stable. There are many equivalent conditions for enforcing
stability of linear systems, including the well-known Lyapunov inequality
$P-{\cal A}^{\top}P{\cal A}\succ 0$ for some $P\succ 0$, and the recently
proposed parameterization in [3], which was used to train stable Koopman
operators for fixed observables in [12]. However, solving optimization
problems with these constraints in an efficient manner is non-trivial,
especially when jointly searching for the observables.

In the following, we present an unconstrained parameterization of ${\cal A}$,
which is a special case of the direct parameterization approach proposed in
[22].

###### Proposition 1

Consider the parametric matrix ${\cal A}(L,R)$ defined as:

${\cal A}(L,R)=2(M_{11}+M_{22}+R-R^{\top})^{-1}M_{21},$ (21)

where

$M:=\begin{bmatrix}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{bmatrix}=LL^{\top}+\epsilon I,$ (22)

with $\epsilon$ a small positive constant. Then for any real-valued
$L\in\mathbb{R}^{2N\times 2N}$ and $R\in\mathbb{R}^{N\times N}$, ${\cal
A}_{0}={\cal A}(L,R)$ is a necessary and sufficient condition for ${\cal
A}_{0}$ to be Schur stable.

###### Proof:

Let $E=(M_{11}+M_{22}+R-R^{\top})/2$, $F=M_{21}$ and $P=M_{22}$. Then we have
${\cal A}(L,R)=E^{-1}F$ and

$M=\begin{bmatrix}E+E^{\top}-P&F^{\top}\\\ F&P\end{bmatrix}.$ (23)

It has been shown that $M\succ 0$ is necessary and sufficient for $E^{-1}F$ to
be Schur stable [26]. Since our parameterization $M=LL^{\top}+\epsilon I$ is
positive definite by construction, this proves sufficiency for ${\cal A}(L,R)$
to be Schur stable. Additionally, all $M$ can be constructed from $L$, e.g.
via Cholesky factorization, and by extension, all $E$, $F$ and $P$ can be
constructed from $L$ and $R$. This completes the proof. ∎

###  IV-C Overall Koopman Model

It may be helpful to think of our model as a linear system with an output
$\hat{x}$ that is an estimate of the original state:

$\displaystyle z(t)$ $\displaystyle={\cal A}z(t-1),$ (24)
$\displaystyle\hat{x}(t)$ $\displaystyle=\phi^{L}(z(t)),$

where $z(0)=\phi(x_{0})$. This system is equivalent to (2). In this form, it
is clear that as long as ${\cal A}$ is stable and $\phi^{L}$ is uniformly
bounded, then the output $\hat{x}$ will always converge to a single
equilibrium.

To summarize, our model parameters consist of:

$\theta=\\{\theta_{NN},\theta_{L},L,R\\}.$ (25)

##  V Learning Framework

###  V-A Optimization Problem

To fit the model parameters (25) to data, we consider the problem of
minimizing the simulation error in the embedding space:

$J_{se}:=\frac{1}{T}\sum_{t=0}^{T}\left\lvert\tilde{z}_{t}-{z}_{t}\right\rvert^{2}_{2},$
(26)

where $\tilde{z_{t}}=\phi(\tilde{x}_{t})$, and $z_{t}={\cal
A}^{t}\phi(\tilde{x}_{0})$. While we could also minimize the simulation error
in $x$, in practice we found this produced poor results.

The complete optimization problem is:

$\min_{\theta\in\Theta}\frac{1}{T}\sum_{t=0}^{T}\left\lvert\phi(\tilde{x}_{t})-{\cal
A}(L,R)^{t}\phi(\tilde{x}_{0})\right\rvert^{2}_{2}+\alpha J_{rec}.$ (27)

The reconstruction loss $J_{rec}$ is defined as

$J_{rec}=\frac{1}{T}\sum_{t=0}^{T}\left\lvert\tilde{x}_{t}-\phi^{L}(\phi(\tilde{x}_{t}))\right\rvert^{2}_{2}.$
(28)

Minimizing $J_{rec}$ gives us an approximate left-inverse $\phi_{L}$ for the
Koopman mapping. The loss $J_{rec}$ can be thought of as a penalty term that
relaxes the constraint

$x=\phi^{L}(\phi(x))\ \forall x,$

and the constant $\alpha$ is a hyperparameter that determines the weighting of
the penalty.

We emphasize two important properties of Problem (27). First, it is an
unconstrained optimization problem. The parameter set $\Theta$ is the space of
real numbers of the appropriate dimensionality. Second, there exists a
differentiable mapping from the parameters $\theta$ to the objective for any
choice of differentiable mapping $\phi_{\theta}$, e.g. using our
parameterization (20) with $\varphi_{\theta}$ as a neural network.

These two properties enable us to find a local optimum to Problem (27) using
any off-the-shelf first-order optimizer in conjunction with an automatic
differentiation (autodiff) toolbox. This significantly simplifies the
implementation of our framework. Using an autodiff software package, one only
needs to write code that evaluates the objective function at each iteration of
the optimization process, and the gradients w.r.t. $\theta$ are automatically
computed via the chain rule. In contrast, constrained problems such as the one
proposed in [12] require specialized algorithms to solve. Although the
objective (27) is nonconvex, deep learning methods have been shown to be
effective at finding approximate global minima for such problems; see [23,
Chapter 21] for example.

###### Remark 2

It is worth noting that our model class is agnostic to the optimization
problem. In fact, the model can be optimized for any differentiable objective
function. This is another advantage of an unconstrained parameterization.

###  V-B Implementation Details

We implemented our learning framework in
PyTorch222https://github.com/pytorch/pytorch and used the Adam optimizer [6]
to solve Problem (27). The neural network parameters $\theta_{NN}$ and
$\theta_{L}$ are initialized using the default scheme in PyTorch, while $L$,
$R$, and $b$ are initialized randomly from a uniform distribution.

####  V-B1 Fast matrix power computation

As explained in Section V-A, the only code we need to implement for solving
Problem (27) is the evaluation of the objective function, which is also the
main computational bottleneck. In particular, repeatedly computing the matrix
power ${\cal A}^{t}$ for the same ${\cal A}$ and many $t$’s can be
computationally inefficient. Here we describe a simple trick to speed up
matrix power computations. Consider the eigendecomposition of ${\cal A}$ given
by $V\Lambda V^{-1}$, where the columns of $V$ are the eigenvectors and
$\Lambda$ is a diagonal matrix of the eigenvalues. Then it is clear that

${\cal A}^{t}=(V\Lambda V^{-1})^{t}=V\Lambda^{t}V^{-1}$ (29)

for integer $t$. Notice that $\Lambda^{t}$ can be computed element-wisely for
each eigenvalue on the diagonal, which offers a significant speed-up over
computing a matrix power. This trick assumes ${\cal A}$ is diagonalizable, but
this can easily be verified in code and, if the condition is not satisfied,
the original matrix power computation can be performed instead.

##  VI Continuous-time Case

In this section, we briefly present the CT formulation of our learning
framework. For an autonomous system governed by an ordinary differential
equation (ODE) $\dot{x}=f(x)$, there exists a semigroup of Koopman operator
$\mathcal{K}^{t}$ associated with the flow map $X(x,t)$ of the system, defined
as:

$\mathcal{K}^{t}\phi(x(t))=\phi(X(x,t)).$ (30)

We refer to the infinitesimal generator of this semigroup as the continuous-
time Koopman operator $\tilde{\mathcal{K}}$ [28]:

$\tilde{\mathcal{K}}\phi(x(t))=\frac{d}{dt}\phi(x(t))=\nabla\phi\cdot
f(x(t)).$ (31)

######  Definition 5 (CT Koopman model)

The continuous-time Koopman model is given by:

$x(t)=\phi^{L}(\exp(A_{\theta}t)\phi(x_{0})),$ (32)

where the Koopman mapping $\phi$ is parameterized as in (20), and the finite-
dimensional matrix $A_{\theta}$ is parameterized as:

$A_{\theta}=(NN^{\top}+\epsilon I)^{-1}(-QQ^{\top}-\epsilon I+{1\over
2}(R-R^{\top})),$ (33)

with parameters $N$, $Q$ and $R$. This is an unconstrained parameterization of
all CT stable (Hurwitz) matrices.

Given full-state trajectory data $\\{\tilde{x}_{k}\\}^{K}_{k=0}$ with
corresponding time $\\{t_{k}\\}^{K}_{k=0}$, we would like to minimize the
simulation error of the Koopman model. The optimization problem is

$\min_{\theta\in\Theta}\frac{1}{K}\sum_{k=0}^{K}\left\lvert\phi_{\theta}(\tilde{x}_{k})-\exp(A_{\theta}t_{k})\phi_{\theta}(\tilde{x}_{0})\right\rvert^{2}_{2}+\alpha
J_{rec},$ (34)

where $J_{rec}$ is as defined in Equation (28). Problem (34) is an
unconstrained optimization problem just like the discrete-time problem, hence
a local minimum can be obtained using a first-order optimizer and an autodiff
software package. A trick similar to that described in Section V-B1 can be
used to compute the matrix exponential in the objective in (34).

Note that in terms of the data required, the only difference between the DT
and CT learning frameworks is that the CT case requires the time corresponding
to each data point. The CT problem can be useful to consider when the data is
sampled at non-uniform time intervals, or when the sampling rate differs
between the training and test scenarios.

##  VII Numerical Examples

TABLE I: Comparison of model sets for our method and prior works. Method |  Learns observables or eigenfunctions |  Continuous or discrete time |  Stability constraint
---|---|---|---
SOC [12] |  Neither |  Discrete |  ✓
LKIS [25] |  Observables |  Discrete |  ✗
[11] |  Eigenfunctions |  Discrete |  ✗
[19] |  Eigenfunctions |  Continuous |  ✓
Ours |  Observables |  Both |  ✓

We validated our framework on the LASA handwriting dataset [5], which consists
of human-drawn trajectories of various letters and
shapes333https://cs.stanford.edu/people/khansari/download.html. It has been
widely used as a benchmark for learning contracting dynamics in continuous-
time [1, 5, 16, 17, 20]. In our results, we trained discrete-time models in
order to compare them with existing DT Koopman learning frameworks.
Contraction is an important constraint for this data set as unconstrained
models can have spurious attractors [5], leading to poor generalization to
unseen initial conditions.

For each shape in the dataset, we attempted to train a discrete-time model
that would regulate to the desired equilibrium point from any initial
condition. To prepare the data for learning DT models, we fitted splines to
the trajectories and re-sampled the datapoints at a uniform time interval. The
state vector was chosen to be
$\tilde{x}_{t}=[y_{t}^{\top},\dot{y}_{t}^{\top}]^{\top}\in\mathbb{R}^{4}$,
where $y_{t}$ and $\dot{y}_{t}$ are the position and velocity vectors at time
$t$. All data was scaled to the range $[-1,1]$ before training. For each shape
in the dataset, we performed leave-one-out cross validation. Test trajectories
are plotted in Figure 3 as solid black lines for a subset of the shapes in the
dataset.

The metric we used to compare different methods was normalized simulation
error (NSE), defined as:

$NSE=\frac{\sum_{t=0}^{T}\left\lvert\hat{x}_{t}-\tilde{x}_{t}\right\rvert^{2}_{2}}{\sum_{t=0}^{T}\left\lvert\tilde{x}_{t}\right\rvert^{2}_{2}},$
(35)

where $\\{\hat{x}\\}^{T}_{t=0}$ is the simulated trajectory using the learned
model, and $\\{\tilde{x}\\}^{T}_{t=0}$ is the true trajectory. The aim of our
comparisons was to evaluate our framework against prior methods for learning
Koopman models. The key differences of some recent frameworks are summarized
in Table I. We did not compare against [11] as they assume some prior
knowledge about the spectrum of the system, which differs from the problem
setting we consider. Due to space constraints, we leave comparisons of CT
frameworks to future work.

In the following, we refer to our framework as SKEL (Stable Koopman Embedding
Learning).

###  VII-A Comparison with other Koopman matrix parameterizations

We compared SKEL against two recently-proposed Koopman learning frameworks,
namely SOC [12] and LKIS [25]. In particular, we compared our unconstrained
stable parameterization of the Koopman operator against a constrained stable
parameterization (SOC), and a unconstrained parameterization without stability
guarantees (LKIS).

The SOC parameterization is given by ${\cal A}=S^{-1}OCS,$ where $O$ is
orthogonal and $C$ is positive-semidefinite with $\lVert C\rVert\leq 1$. A
projected gradient descent method was used to solve the optimization problem.
The LKIS parameterization is ${\cal A}=Y_{1}Y_{2}^{\dagger},$ where $Y_{1}$
and $Y_{2}$ are as defined in Equation (6), with parametric $\phi(x)$.

To make it a fair comparison, we kept all other aspects of the optimization
problem the same, i.e. using simulation error as the optimization objective
and using parametric observables of the form (20). We were interested mainly
in comparing parameterizations of the Koopman operator as our framework is
agnostic to choice of objective and observables, and these choices often
depend on the particular application.

All instances of $\varphi(x)$ were fully-connected feedforward neural networks
with ReLU (rectified linear units) activation functions, 2 hidden layers with
50 nodes each and an output dimensionality of 20. Hyperparameter values were
chosen to be $\alpha=10^{3}$ and $\epsilon=10^{-8}$.

A boxplot of the normalized simulation error for the three methods is shown in
Figure 1. It is clear that SKEL achieves the lowest median NSE on the test set
with 95% confidence. From Figure 2, it can be seen that LKIS actually attains
the lowest training error, but does not generalize to the test set as well as
SKEL. This can be seen as a symptom of overfitting, and shows that the
stability guarantees of SKEL have a regularizing effect on the model. With
regards to SOC, we observed that the constrained optimization problem would
often converge to poor local minima, which is reflected in the relatively high
training and test errors.

Figure 1: Comparison of SKEL with other Koopman learning methods. Outliers
were clipped for better visibility of boxes. Number of outliers with NSE $>1$
from left to right: 1 (SKEL), 15 (LKIS), 0 (SOC). Figure 2: Training loss (Eq.
(27)) for each method

(a) SKEL (ours)

(b) LKIS [25]

Figure 3: Simulations of SKEL and LKIS models on test data. Trajectories from
the models are shown as red dotted lines, while the true trajectory is shown
as a solid black line. Initial conditions were sampled from a square region of
width 2mm centered at the start point of the true trajectory. The target point
is marked by a black star.

###  VII-B Robustness to perturbations in initial condition

We performed a qualitative evaluation of the robustness of the models to small
perturbations in the initial condition of the test trajectory. We compared
only SKEL and LKIS as it was clear from Figure 1 that SOC underperformed in
this setting. The results are plotted in Figure 3. It can be seen that the
SKEL models produce trajectories that converge to each other due to their
contracting property, whereas the LKIS models behave unpredictably, indicating
instability of the learned model.

##  VIII Conclusion

We have presented a novel Koopman learning framework that jointly models the
Koopman operator and observables while guaranteeing model stability, via an
unconstrained optimization problem. We showed that our framework outperforms
existing Koopman methods on a real-world handwriting problem and achieves the
lowest median simulation error. Further work can be done to extend this
framework to controlled systems.

## References

  * [1] C. Blocher, M. Saveriano, and D. Lee.  Learning stable dynamical systems using contraction theory.  In 14th Int. Conf. on Ubiquitous Robots and Ambient Intell. (URAI), pages 124–129, 2017.
  * [2] L. Brivadis, V. Andrieu, and U. Serres.  Luenberger observers for discrete-time nonlinear systems.  In the 58th IEEE Conf. on Decision and Control (CDC), pages 3435–3440. IEEE, 2019.
  * [3] N. Gillis, M. Karow, and P. Sharma.  A note on approximating the nearest stable discrete-time descriptor systems with fixed rank.  Appl. Numerical Math., 148:131–139, 2020.
  * [4] M. Haseli and J. Cortes.  Learning koopman eigenfunctions and invariant subspaces from data: Symmetric subspace decomposition.  IEEE Trans. on Automatic Control, pages 1–1, 2021.
  * [5] S. M. Khansari-Zadeh and A. Billard.  Learning stable nonlinear dynamical systems with gaussian mixture models.  IEEE Trans. on Robotics, 27(5):943–957, 2011.
  * [6] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [7] B. O. Koopman.  Hamiltonian systems and transformation in Hilbert space.  Proc. of the National Academy of Sciences, 17(5):315–318, May 1931\.
  * [8] M. Korda and I. Mezić.  Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control.  Automatica, 93:149 – 160, 2018.
  * [9] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis.  Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator.  Chaos: An Interdisciplinary J. of Nonlinear Science, 27(10):103111, 2017.
  * [10] W. Lohmiller and J.-J. E. Slotine.  On contraction analysis for non-linear systems.  Automatica, 34(6):683–696, 1998.
  * [11] B. Lusch, J. N. Kutz, and S. L. Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature Comm., 9(1):4950, Nov. 2018.
  * [12] G. Mamakoukas, O. Xherija, and T. Murphey.  Memory-efficient learning of stable linear dynamical systems for prediction and control.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Syst., volume 33, pages 13527–13538. Curran Associates, Inc., 2020.
  * [13] I. R. Manchester, M. Revay, and R. Wang.  Contraction-based methods for stable identification and robust machine learning: a tutorial.  In the 60th IEEE Conf. on Decision and Control (CDC), 2021.
  * [14] A. Mardt, L. Pasquali, H. Wu, and F. Noé.  Vampnets for deep learning of molecular kinetics.  Nature Comm., 9(1):5, Jan. 2018.
  * [15] A. Mauroy and J. Goncalves.  Koopman-based lifting techniques for nonlinear systems identification.  IEEE Trans. on Automatic Control, 65(6):2550–2565, 2020.
  * [16] S. Mohammad Khansari-Zadeh and A. Billard.  Learning control lyapunov function to ensure stability of dynamical system-based robot reaching motions.  Robotics and Autonomous Syst., 62(6):752–765, 2014.
  * [17] K. Neumann, A. Lemme, and J. J. Steil.  Neural learning of stable dynamical systems based on data-driven lyapunov candidates.  In IEEE/RSJ Int. Conf. on Intell. Robots and Syst., pages 1216–1222, 2013.
  * [18] S. E. Otto and C. W. Rowley.  Linearly recurrent autoencoder networks for learning dynamics.  SIAM J. on Applied Dynamical Syst., 18(1):558–593, 2019.
  * [19] S. Pan and K. Duraisamy.  Physics-informed probabilistic learning of linear embeddings of nonlinear dynamics with guaranteed stability.  SIAM J. on Appl. Dynamical Syst., 19(1):480–509, 2020.
  * [20] H. Ravichandar, I. Salehi, and A. Dani.  Learning partially contracting dynamical systems from demonstrations.  In S. Levine, V. Vanhoucke, and K. Goldberg, editors, Proc. of the 1st Annual Conf. on Robot Learning, volume 78 of Proceedings of Machine Learning Research, pages 369–378. PMLR, 13–15 Nov 2017.
  * [21] M. Revay, R. Wang, and I. R. Manchester.  A convex parameterization of robust recurrent neural networks.  IEEE Control Syst. Lett., 5(4):1363–1368, 2021.
  * [22] M. Revay, R. Wang, and I. R. Manchester.  Recurrent equilibrium networks: Unconstrained learning of stable and robust dynamical models.  2021\.  arXiv:2104.05942.
  * [23] T. Roughgarden.  Beyond the Worst-case Analysis of Algorithms.  Cambridge University Press, 2020.
  * [24] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  J. of Fluid Mechanics, 656:5–28, 2010.
  * [25] N. Takeishi, Y. Kawahara, and T. Yairi.  Learning koopman invariant subspaces for dynamic mode decomposition.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 1130–1140. Curran Associates, Inc., 2017.
  * [26] M. M. Tobenkin, I. R. Manchester, and A. Megretski.  Convex parameterizations and fidelity bounds for nonlinear identification and reduced-order modelling.  IEEE Trans. on Automatic Control, 62(7):3679–3686, 2017.
  * [27] J. Umenberger and I. R. Manchester.  Convex bounds for equation error in stable nonlinear identification.  IEEE Control Syst. Lett., 3(1):73–78, 2019.
  * [28] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley.  A data–driven approximation of the koopman operator: Extending dynamic mode decomposition.  J. of Nonlinear Science, 25(6):1307–1346, Dec 2015.
  * [29] E. Yeung, S. Kundu, and N. Hodas.  Learning deep neural network representations for koopman operators of nonlinear dynamical systems.  In American Control Conf. (ACC), pages 4832–4839, 2019.
  * [30] B. Yi and I. R. Manchester.  On the equivalence of contraction and Koopman approaches for nonlinear stability and control.  In the 60th IEEE Conf. on Decision and Control (CDC), 2021.  (See ArXiv:2103.15033 for an extended version).


# Numerical studies of triangulated vesicles with anisotropic membrane
inclusions

Yoav Ravid, Nir Gov
Department of Chemical and Biological Physics
Weizmann Institute of Science
Rehovot 7610001, Israel

&Mitja Drab, Aleš Iglič
Laboratory of Physics
Faculty of Electrical Engineering
University of Ljubljana, 1000 Ljubljana, Slovenia
Corresponding author<EMAIL_ADDRESS>
Veronika Kralj-Iglič
Laboratory of Clinical Biophysics
Faculty of Health Sciences
University of Ljubljana, 1000 Ljubljana, Slovenia
&Samo Penič
Laboratory of Bioelectromagnetics
Faculty of Electrical Engineering
University of Ljubljana, 1000 Ljubljana, Slovenia

###### Abstract

In this study, we implement the deviatoric curvature model to examine
dynamically triangulated surfaces with anisotropic membrane inclusions. The
Monte-Carlo numerical scheme is devised to not only minimize the total bending
energy of the membrane but also the in-plane nematic order of the inclusions
by considering the mismatch between the curvature of the membrane and the
intrinsic curvature of the inclusion. Neighboring inclusions can either
attract with nearest-neighbor interaction or with a nematic interaction
derived from liquid crystal theory. Orientational order determines whether
vesicles fully covered with inclusions result in bulbs connected by necks or
long tubes. Remarkably, when inclusions on vesicles with no vacancies interact
non-nematically, a spontaneous local order can lead to a bulb transition which
may have implications in cell or organelle division. Furthermore we find that
average nematic order is inversely proportional to the number of thin necks
formed in the vesicles. Our method shows good convergence and is suitable for
further upgrades, for example to vesicles constrained by volume.

_K_ eywords Anisotropic $\cdot$ Membrane inclusions $\cdot$ Cell division
$\cdot$ Monte-Carlo $\cdot$ Vesicles

##  1 Introduction

Cellular shape holds immense importance in various functions like division,
movement, and signaling. In biophysics, there’s a significant interest in
comprehending how cells take shape. Research involving theories and computer
simulations indicates that the physical traits of the cell membrane, like its
uneven components and dynamic forces, significantly influence structure. The
precise interplay between these elements in shaping cells remains uncertain.

Membrane shape depends on the intrinsic shape of the membrane’s molecular
constituents [1, 2, 3, 4, 5] and their interactions with other components,
such as membrane skeleton and cytoskeleton [6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16]. It has been shown that a non-homogeneous lateral distribution and
phase separation of membrane inclusions may be a driving force of cell shape
transformations and necessary for the stabilization of highly curved membrane
structures [17, 18, 19, 20, 21, 22, 23, 2, 3].

In recent years, there have been significant advances in experimental
techniques for studying cell mechanics [24, 25, 26]. However, theoretical and
computational studies remain important for gaining insight into the underlying
physical mechanisms. For example, the study of vesicle shapes driven by
coupling curved inclusions and active cytoskeletal forces has shown that the
inhomogeneous nature of membrane constituents, such as curved inclusions, can
have a significant impact on vesicle shapes [11, 20, 27, 28, 29, 30, 31, 32,
33].

In this work, we apply the deviatoric curvature model to investigate
dynamically triangulated surfaces containing anisotropic membrane inclusions
that are banana shaped. A Monte-Carlo numerical approach is employed not only
to reduce the overall bending energy of the membrane but also to minimize the
in-plane nematic order of the inclusions. This is achieved by considering the
discrepancy between the curvature of the membrane and the inclusion it
contains. We find that the orientational order has an effect on the overall
limiting cell shape. Additionally, we present a possible mechanism of cell
division when inclusions interact non-nematically through spontaneous
ordering. Furthermore, we find that the total bonding energy of the fully
covered vesicles is inversely proportional to the number of necks and buds in
the vesicles.

##  2 Theoretical background

Theoretical background is based on the derivation of membrane energy and
anisotropic inclusion interactions, as outlined by Kralj-Iglič et al. [34, 2,
4]. The model presents the membrane as a 2D anisotropic continuum surface. The
elastic energy of a small membrane section becomes null when its principal
curvatures $C_{1}$ and $C_{2}$ align with their intrinsic counterparts of the
anisotropic inclusions $C_{1m}$ and $C_{2m}$, and when the orientations of the
local membrane curvature tensor $C$ coincide with the intrinsic membrane
curvature tensor $C_{m}$. A shape with these properties throughout exhibits
zero elastic energy.

In essence, the elastic energy per unit area of a thin plate element with area
$dA$ is defined by the discrepancy between its actual local membrane curvature
and its intrinsic curvature. Representing the actual and intrinsic curvatures
with tensors $C$ and $C_{m}$, respectively, where in their principal systems,
the curvature tensor matrices feature only diagonal elements [2, 34]:

$C=\begin{pmatrix}C_{1}&0\\\ 0&C_{2}\end{pmatrix},\quad
C_{m}=\begin{pmatrix}C_{1m}&0\\\ 0&C_{2m}\end{pmatrix}.$

In this examination, the principal systems of tensors generally undergo an
angle $\omega$ rotation in the tangent plane. The disparity between the actual
and intrinsic local membrane curvature is expressed by the tensor [2, 4]
$M=RC_{m}R^{-1}-C$, involving the rotation matrix

$R=\begin{pmatrix}\cos\omega&-\sin\omega\\\
\sin\omega&\cos\omega\end{pmatrix}.$

A small shell patch adjusts to conform to the actual membrane, reflecting the
energy required for deformation (see Figure 1(b)). The approximate elastic
energy $E_{1}$ unfolds through an expansion in powers of all independent
invariants of the tensor $M$. Employing the trace and determinant of the
tensor as invariants yields [34, 2, 4, 3]:

$E_{1}=\int\frac{K_{1}}{2}(\mathrm{Tr}M)^{2}+K_{2}\mathrm{Det}M\,dA,$ (1)

integrated across the entire vesicle. With constants $K_{1}$ and $K_{2}$, the
energy $E_{1}$ is [34, 2, 4, 35]:

$E_{1}=\int(2K_{1}+K_{2})(H-H_{m})^{2}-K_{2}\left(D^{2}-2DD_{m}\cos
2\omega+D_{m}^{2}\right)\,dA,$ (2)

where $D=(C_{1}-C_{2})/2$ stands as the curvature deviator, an invariant of
the curvature tensor
($D^{2}=(\mathrm{Tr}(C)/2)^{2}-\mathrm{Det}(C)=H^{2}-C_{1}C_{2}$);
$H_{m}=(C_{1m}+C_{2m})/2$ denotes the intrinsic mean curvature; and
$D_{m}=(C_{1m}-C_{2m})/2$ represents the intrinsic curvature deviator.

In the case of an isotropic membrane where $D_{m}=0$, it becomes evident that
this equates to the Helfrich bending energy density [36] described by
$E_{b}=k^{2}_{c}(2H-C_{0})^{2}+k_{G}K$, with $H=(C_{1}+C_{2})/2$ as the mean
curvature, $C_{0}$ as the spontaneous curvature, $K=C_{1}C_{2}$ as the
Gaussian curvature, and $k_{c}$ and $k_{G}$ as the bending and splay modulus,
respectively [34, 3, 4].

Figure 1: Example of a banana inclusion. For this particular inclusion,
$H_{m}=D_{m}=0.25$ (a). The mismatch between a membrane inclusion and local
curvature of the triangulated membrane (b). Here, $\vec{n}$ is the orientation
of the inclusion and $\vec{t}$ is its perpendicular direction in the local
tangent plane.

###  2.1 Binding energy

Inclusions may bond on the membrane by a simple energy gain for two
neighboring vertices.

$E_{2}=\begin{cases}-w&i,j\text{ are binding}\\\ 0&\text{else}\end{cases}$ (3)

where the type and the binding strength $w$ are vertex properties, reflecting
aggregation. Throughout this paper, we refer to this type of bonding as normal
bonding. This kind of bonding is observed in surfactants and proteins, since
it arises from spontaneous self-assembly on the membrane [37, 11].

###  2.2 Nematic energy

If we want to impose a directional energy penalty to neighboring inclusions,
we use a nematic-nematic interaction employed from liquid crystal theory. The
energy between neighboring inclusions will be lower if their principal
curvatures are aligned. We model this based on Frank’s free energy density for
nematic liquid crystals [38]:

$\displaystyle E_{2}$
$\displaystyle=\int\left(\frac{k_{G}}{2}\left(\nabla\cdot\vec{n}\right)^{2}+\frac{k_{c}}{2}\left(\nabla\cdot\vec{t}\right)^{2}\right)r_{v}\,dA$
(4)

Here, $\nabla$ represents the covariant derivative on the curved surface,
$\vec{n}$ denotes the orientation of the inclusion, and $\vec{t}$ signifies
its perpendicular direction in the local tangent plane. The constants $k_{G}$
and $k_{c}$ correspond to the splay and bending elastic constants governing
in-plane nematic interactions. The variable $r_{v}$ takes a value of 1 if a
nematic is present on a vertex; otherwise, it is 0.

A discrete form of this energy is employed to facilitate implementation in
Monte-Carlo (MC) simulations. If we assume a one-constant approximation
($k_{G}=k_{c}$), equation 4 can be rewritten in a way that makes the
implementation suitable for MC simulations (Lebwohl-Lasher model) [39]:

$E_{2}\approx-\int\sum_{k=1}^{N}\epsilon_{LL}^{k,k}\sum_{i>j}\left(\dfrac{3}{2}(\vec{n_{i}}\cdot\vec{n_{j}})^{2}-\dfrac{1}{2}\right)r_{v}\,dA.$
(5)

Here, $\epsilon_{LL}$ is the strength of the nematic interaction. The sum
$\sum_{i>j}$ is over all the nearest neighbor (i, j) vertices on the
triangulated grid, promoting alignment among the neighbouring orientation
vectors. An even simpler form of this approximation for the in-plane
orientational field is the XY model on a random surface [40]:

$E_{2}\approx-w\int\sum_{i>j}(\vec{n}_{i}\cdot\vec{n}_{j})^{2}\,dA.$ (6)

Here, $w$ represents the strength of the direct interaction constant, and the
summation runs over all inclusion-inclusion pairs. The sum of $E_{1}$ and
$E_{2}$ is minimized numerically, while $E_{2}$ is either normal (equation 3)
or nematic bonding (equation 6). There is no volume restricting the vesicle
shapes. All the inclusions in this work are banana-shaped (see Figure 1.1(a)).
Normal interaction reflects the tendency of inclusions to self-aggregate.
Nematic interaction goes further and implies the tendency to both self-
aggregate and also align inclusions’ principal curvatures.

The details of the Monte-Carlo procedure are given in Appendix 6.1 and the
details of the mesh generation and finding the principal curvatures are given
in in Appendix 6.2.1.

##  3 Results and discussion

###  3.1 Small meshes

We start by examining smaller vesicles ($N=502$ vertices) with a fraction of
their surface $\rho$ covered with banana shaped inclusions as shown in Figure
1(a). One of their main curvatures is positive and another is zero, which
translates to $H_{m}=D_{m}=0.25$. With increasing $\rho$ and normal
interaction between inclusions, we get the sequence shown in Figure 2. Due to
the inclusions’ strong affinity to aggregate ($w=3$), they form a rim around
the perimeter of the vesicle even at low values of $\rho$. With increasing
$\rho$, the vesicles become flattened and oblate and later elongated (but not
prolate). Only above $\rho=0.7$ does the high inclusion content clip the
vesicle into two bulbs connected by a thin neck.

Figure 2: A sequence of equilibrium shapes for a small mesh with normal
interaction between banana inclusions. The constant of inclusion interaction
is $w=3$, while the curvatures are $H_{m}=D_{m}=0.25$. The inclusions are
marked in red, the empty membrane in blue. The principal direction of the
curvature $C_{1m}$ is shown with the white lines. The average order parameter
of each shape $\langle S\rangle$ is given in the bottom and calculated by
equation 8

We want to quantify the degree of order between neighboring anisotropic
inclusions. For this reason we can define the nematic order of the inclusion
at vertex $i$ by

$S_{i}=\frac{1}{2}(3\cos^{2}{\theta}-1),$ (7)

where the angle $\theta$ measures the angle between directors (or principal
curvature $C_{1m}$) of neighboring inclusions on the membrane. The more
aligned the neighboring inclusions are, the more this value approaches unity.
In our numerical approach it is calculated as follows: at vertex $i$, we
average the dot product of the director $\vec{d}_{i}$ with all immediate
neighboring vertices (see Appendix 6.2.3). The average nematic order across
the entire membrane is then

$\langle S\rangle=\frac{1}{N}\sum_{i=1}^{N}S_{i}.$ (8)

Returning to Figure 2, we notice that predictably the nematic order increases
with increasing $\rho$, but peaks just before the bulb transition. Their
binding promotes a rim on the membrane that accomodates their spontaneous
curvature best when they all align in a similar direction (see Figure 2 for
$\rho=$ 0.3–0.7). When the inclusion density is large enough, the ordering
ceases, resulting in two prolate bulbs connected by a neck where nematic order
decreases. We call this the bulb transition. The average order parameter at
the bottom of Figure 1.2. shows that the bulb transition decreases the average
order.

Figure 3: Simulations of vesicles fully covered in banana inclusions
($H_{m}=D_{m}=0.25$) for a small mesh (N=502) and normal interaction
($w=3$)(A) and nematic interaction (B). With consequent MC steps, the shapes
relax to an energy minimum. For normal interactions, the resulting shape is a
two bulbs connected by a thin neck. For the nematic case, the limiting shape
is an elongated and thin cylinder. The average order parameter $\langle
S\rangle$ for both cases is shown below the main graph.

Next we investigate vesicles with all places on the membrane occupied by
inclusions. These two limiting cases of $\rho=1$ are shown in Figure 3. Here,
we observe the evolution of the vesicle shapes from its spherical start. We
may think of consequential MC steps as simulation time shown on the $x$-axis.

First, we look at the normal interaction between inclusions shown in Figure
3(A). The inclusions predictably show no sign of nematic order on the
membrane, prompting neighbours to point in a random direction and hence
promote the spherical shape that persists for a long time. This is because the
average curvature of many random directions results in a net zero curvature
and a spherical shape. However, the fluctuations with time lead again to a
bulb transition. The bulb transition coincides with a sharp decline in the
bending energy (accompanied by a decline in volume) shown in Figure 3(A). The
average nematic order slightly increases before it returns to its pre-
transition values.

Upon closer inspection of an increase in $\langle S\rangle$ preceding the bulb
transition for Figure 3(A), we notice an interesting sequence of events shown
in Figure 4. Inclusions spontaneously order along the circumference of the
prolate vesicle, resulting in an increase of $\langle S\rangle$ (see Figure
4(B)). This is followed by a sudden narrowing of the midsection and results in
a bulb transition (Figure 4(C)). Here, $\langle S\rangle$ and bending energy
decrease relative to their pre-transition values. The spontaneous ordering of
the inclusions can therefore have lasting impact on the morphology of the
vesicle and may facilitate in the mechanism of cell division.

Figure 4: A close-up of the average order parameter during the bulb transition
(see Figure 3(A)). A spontaneous ordering of the inclusions results in the
transition from a sphere-like shape (A) to bulbs connected by a thin neck (B).
After the transition, $\langle S\rangle$ decreases together with the bending
energy.

If we run the same simulation with nematic interactions between inclusions,
the limiting shape is an elongated cylinder shown in Figure 3(B). The bending
energy continuously decreases in proportion with the average radius of the
cylinder and leveling out at a limiting value. The overall bending energy of
the thin cylinder is lower than the bulbs connected by a thin neck.

The nematic order interestingly slightly decreases (cca. 10%) as the cylinder
becomes elongated. This decrease may be a result of inclusions that are packed
around a tighter radius. Let us envision this with a graphical example. In the
limit of a very thin cylinder there may be only 5 vertices encircling the
smaller cylinder as shown in Figure 5. Even though the inclusions’ directions
are aligned, neighboring angles are relatively large resulting in near zero
nematic order (see Figure 5).

Figure 5: Low nematic order as a consequence of very thin necks of meshes.
Although the curvature of the inclusions fits well around the circumference,
the calculation of average nematic order given in equation 7 will be low
because of the large angles ($\theta$) between neighboring inclusions. This
artefact also accounts for near zero nematic order in the necks of buds.

###  3.2 Larger meshes

Now we turn our attention to larger vesicles ($N=2002$). We start again by
gradually increasing the fraction of inclusions in the membrane with other
parameters staying the same as before. The results are given in Figure 6.

Figure 6: Simulation results for larger meshes ($N=2002$), $H_{m}=D_{m}=0.25$,
normal interaction ($w=3$) after $2\cdot 10^{6}$ MC steps in dependence of
inclusion coverage fraction $\rho$. The color code is the same as in Figure 2.

At $\rho=0.2$, the equilibrium shapes resemble quasi-spherical vesicles with
isolated patches of inclusions forming ridges along the membrane surface. The
ridges are formed due to the accommodation of the principal curvature of the
inclusions $C_{1m}$. Increasing $\rho$ results in patches joining together to
form a single ridge around the outer circumference of a flattened vesicle. As
we can see in Figure 6 for $\rho=0.4$, the continuous ridge forces the vesicle
shape away from the oblate pancake into a global saddle shape. Increasing
$\rho$ further results in separation of a single saddle pancake to two
separate ones. With $\rho=0.6$, each of the pancakes starts forming buds
connected by a very thin neck. At $\rho=0.8$, nearly the whole surface of the
mesh consists of buds covered with inclusions to resemble a pearling shape.

Figure 7: Nematic order as a heat map for small (A, B) and larger meshes (C,
D). When the inclusion interaction is normal (A, C), there is little neighbor
ordering shown by values close to 0 and a blue color. Conversely, red color
represents strong nematic order as expected by the nematic neighbor
interaction (B, D). Some spontaneous ordering may happen as seen by the orange
patch in (A). The heatmap $S_{i}$ gives the local order parameter.

Figure 7 shows a heat map of the nematic order for a fully covered vesicle for
a small (A, B) and larger mesh (C, D). Normal and nematic interaction between
inclusions is considered for (A, C) and (B, D), respectively. We observe that
nematic ordering results in elongated cylinders with nematic order near unity
everywhere except at the poles of the cylinders, where a nematic defect is
observed in blue (Figure 7(B)). For larger meshes with nematic interactions
the budding from the main trunk of the vesicle still occur, but the buds
themselves are cylindrical. As expected, we notice that the nematic order is
close to 0 in the necks where the radius is very thin.

Figure 8:

Figure 9:

Figure 10:

Figure 11:

For each final shape at constant $\rho$ for Figures 6, Figures 11–11 show the
bending energy, bonding energy and volume as a function of MC steps of the
simulation. We observe that the vesicles approach their final volume faster in
proportion to $\rho$. This likely happens as a consequence of budding which is
more prevalent at higher inclusion content. The limiting pearling shape has
the least volume per given vesicle area and the largest bending energy of all
the shapes in the sequence of Figure 6.

###  3.3 Low nematic order is porportional to the number of buds

Lastly we look at fully covered vesicles with nematic interaction with
deviatoric curvature $H_{m}=D_{m}=0.75$ and all other parameters same as
before. The banana shape of the inclusion is preserved, but its curvature is
increased. Figure 12(A) shows the graphs of bending and bonding energy with
respect to $w$.

Figure 12: Total energy per vertex in dependence on $w$ for $H_{m}=D_{m}=0.75$
(A). Highly curved anisotropic inclusions result in very thin spaghetti-like
shapes for normal interaction ($w=2$). The nematic order here is low because
of the thin radius of the tubes (B).

Bonding energy depends on two things; the strength of interaction $w$ and
local nematic order. We find that high bonding energy (and hence worst nematic
order) goes to membranes with highest curved inclusions $H_{m}=0.75$. This
implies at least one of two things are taking place: either vesicles form many
offshoot buds with their respective necks, or the vesicles resemble very
elongated thin tubular structures. Figure 12(B) shows such a final shape
($H_{m}=D_{m}=0.75$) and confirms our prediction of a low nematic order (and
high bonding energy).

##  4 Conclusion

Membrane deformations in cells occur due to curved inclusions working in
synergy with each other. Cell mechanics studied from a simple model with few
assumtions might shed new light on formation of vesicle morphology.

In this paper we demonstrated a new way of simulating anisotropic membrane
inclusions by implementing the deviatoric elasticity model that was developed
nearly twenty years ago [34]. Although our investigation was limited only to
banana membrane inclusions, we found that orientational order determines
whether vesicles fully covered with inclusions result in pearls connected by
necks or long tubes with few necks. Additionally, we observed that the average
nematic order is inversely proportional to the number of thin necks formed by
the vesicles.

The dynamics shown in Figure 4 may indicate a mechanism for primitive cell
division or division of organelles inside the cell, using many curved banana-
shaped membrane proteins. It is crucial that the membrane inclusions interact
non-nematically. Nematic ordering of inclusions result in long narrow tubes
and does not lead to fission-like geometries. This is a mechanism that will be
explored in further papers.

An obvious advantage of this model is the absence of limitations of axial
symmetry. On the other hand, the model has few parameters, and its phase space
offers further investigation. One could for example look at inclusion
concentrations lower than $10\%$ or use a combination of curved inclusions.
Some processes, like clathrin mediated endocytosis and formation of three way
junctions in the endoplasmic reticulum [41, 42], are examples of membrane
deformation and stabilisation of induced curvature due to a combination of
curved inclusions.

A further development would entail a constraint of constant volume. This could
be achieved by introducing a new term to the Hamiltonian to include a volume
potential to force the vesicles toward a wanted set volume. In this way,
vesicle morphologies could be compared to the standard empty Helfrich diagram
[36] and study the effects that inclusions have on the overall morphology.
Furthermore, membranes with active inclusions could be modelled by adding a
non-conservative force acting perpendicular away from the vesicle to simulate
the polymerization of actin [11, 43].

##  5 Acknowledgements

This work was supported by the Slovenian Research Agency (ARIS) through Grants
Nos. J3-3066 and J2-4447 and Programme No. P2-0232. N. G. is the incumbent of
the Lee and William Abramowitz Professorial Chair of Biophysics, and
acknowledges support by the Ben May Center for Theoryand Computation, and the
Israel Science Foundation (Grant No. 207/22).

##  6 Appendix - Methods

###  6.1 Monte-Carlo procedure

The membrane is represented by a set of $N$ vertices that are linked by
tethers of variable length $l$ to form a closed, dynamically triangulated,
self-avoiding two-dimensional network of approximately 2N triangles and with
the topology of a sphere [44, 45]. The lengths of the tethers can vary between
a minimal and a maximal value, $l_{min}$, and $l_{max}$, respectively. Self-
avoidance of the network is ensured by choosing the appropriate values for
$l_{max}$ and the maximal displacement of the vertex $s$ in a single updating
step.

One Monte-Carlo sweep (MCs) consists of individual attempts to displace each
of the $N$ vertices by a random increment in the sphere with radius $s$,
centered at the vertex, followed by $R_{B}N$ attempts to flip a randomly
chosen bond. We denote $R_{B}$ as the bond-flip ratio, which defines how many
attempts to flip a bond are made per one attempt to move a vertex in one MCs.
Note that the bond-flip ratio is connected to the lateral diffusion
coefficient within the membrane, i.e. to the membrane viscosity. In this work
we have chosen $R_{B}=3$, $s/l_{min}=0.15$ and $l_{max}/l_{min}=1.7$. The
dynamically triangulated network acquires its lateral fluidity from a bond
flip mechanism. A single bond-flip involves the four vertices of two
neighboring triangles. The tether connecting the two vertices in diagonal
direction is cut and reestablished between the other two, previously
unconnected, vertices. The self-avoidance of the network is implemented by
ensuring that no vertex can penetrate through the triangular network and that
no bond can cut through another bond [46, 47].

###  6.2 Anisotropic code details

####  6.2.1 Representing the membrane as a mesh

Our solver is called Trisurf. It models the vesicle as a closed, triangulated
surface: a graph with vertices $i\in V$ and edges $e_{ij}\in E$, and an
auxiliary set of triangles $t_{ijk}\in T$. The triangles make the
approximation of the surface, but it is the vertices which are the principle
dynamical entities which holds the properties of the membrane (intrinsic
curvature, membrane composition, nematic director, etc.) and move in space.

From the position of the vertices, $\vec{x}_{i}$ we can compute the normal
vector to each triangle $\vec{N}_{ijk}$ and circumcenter $\vec{O}_{ijk}$
(center of the circle containing the position of
$\vec{x}_{i},\vec{x}_{j},\vec{x}_{k}$). This allows us to divide the triangle
to six parts and assign two pieces to each vertex.

For vertex $i$, these are the sub-triangle between the vertex position
$\vec{x}_{i}$, the triangle center $\vec{O}_{ijk}$, and the middle of one of
the edges $\frac{\vec{x}_{i}+\vec{x}_{j}}{2}$, and a similar sub-triangle for
the other edge $ik$. We can denote the vector for the side between the vertex
and the edge middle as half the edge length
$\vec{e}_{ij}=\vec{x}_{j}-\vec{x}_{i}$ and the side between center and the
edge middle, which is half of the dual (voronoi) edge, as
$\vec{\sigma}^{i}_{jk}$ (the other half of the voronoi edge is on the
neighboring triangle $\sigma^{i}_{\ell j}$). We can compute this $\sigma$ from
the circumcenter and the middle of $\vec{x}_{i},\vec{x}_{j}$:

$\vec{\sigma}^{i}_{jk}=\frac{\vec{x}_{i}+\vec{x}_{j}}{2}-\vec{O}_{ijk}$ (9)

With this, we assign an area $A\left(i\right)$ and a normal $N\left(i\right)$
for each vertex $i$, by running over the neighbors $j-1,j,j+1\ldots$ and the
adjacent triangles $\left(i,j-1,j\right),\left(i,j,j+1\right)\ldots$

$A\left(i\right)=\sum_{\left\langle
i,j\right\rangle}\frac{1}{2}\frac{\left|e_{ij}\right|}{2}\left|\sigma^{i}_{j,j+1}\right|+\frac{1}{2}\frac{\left|e_{ij}\right|}{2}\left|\sigma^{i}_{j,j-1}\right|$
(10) $\vec{N}\left(i\right)=\frac{\sum_{\left\langle
i,j\right\rangle}\vec{N}_{ij,j+1}\left|e_{ij}\right|\left|\sigma^{i}_{j,j+1}\right|+\vec{N}_{i,j-1,j}\left|e_{ij}\right|\left|\sigma^{i}_{j,j-1}\right|}{\left|\ldots\right|}$
(11)

Where the $j$ vertices are the counterclockwise ordered neighbors of $i$.

The fluidity of the surface is achieved by bond-flips, where a bond $ij$ and
the two triangles that share it is $ijk,ji\ell$ are replaced by a cross bond
$k\ell$ and two triangles $i\ell k,jk\ell$, which allows vertices to change
neighbors.

####  6.2.2 Anisotropic curvature on the vesicle

To get the anisotropic bending energy of the surface, we use the method by
[40] to estimate the shape operator matrix $S$ on each vertex, which
represents the shape of the surface at the point. We then calculate the
mismatch tensor $M=S-C_{m}$, where $C_{m}$ is the intrinsic curvature tensor
whose direction is determined by the director and the other tangent vector to
the local normal ($\hat{t}=\hat{N}\times\hat{d}$)

$C_{m}=\frac{C_{0}+D_{0}}{2}\hat{d}\otimes\hat{d}+\frac{C_{0}-D_{0}}{2}\hat{t}\otimes\hat{t}$
(12)

where $C_{0},D_{0}$ are the spontaneous curvature and spontaneous deviator at
the vertex, respectively, which reflects the physical characteristics of local
membrane composition.

The bending energy is calculated by inserting the mismatch tensor in the
Hamiltonian

$E_{1}=\frac{K_{1}}{2}\left(\mathrm{Tr}M\right)^{2}+K_{2}\mathrm{Det}M$ (13)

Where $K_{1}$ and $K_{2}$ are the bending moduli of the vertex, again
reflecting physical parameters due to local composition.

To calculate the shape curvature of a vertex $i$ based on [40], each edge $ij$
is assigned a shape tensor estimation

$S_{ij}=h_{ij}\;\vec{b}\otimes\vec{b}$ (14)

$\vec{b}$ is the binormal at the edge $\hat{N}_{ij}\times\vec{e}_{ij}$, where
$\vec{e}_{ij}=\vec{x}_{j}-\vec{x}_{i}$ is the edge vector and $\hat{N}_{ij}$
is the normal of the edge
$\hat{N}_{ij}=\left(\hat{N}_{i,j-1,j}+\hat{N}_{i,j,j+1}\right)/\left|\ldots\right|$
which is the sum of the normal of the two triangles sharing the edge,
normalized. $h_{ij}$ is a factor representing the directional derivative of
the area $\approx\nabla_{p}A$

$h_{ij}=2\left|e_{ij}\right|\cos\left(\frac{\Phi}{2}\right)$ (15)

Where $\Phi$ is the dihedral angle (angle between the two triangle sharing the
edge). Luckily there is a simple triple product formula for this factor

$h_{ij}=2\left|\vec{e}_{ij}\right|\vec{N}_{ij}\cdot\left(\vec{N}_{i,j-1,j}\times\hat{e}_{ij}\right)$
(16)

The cross product of the edge direction and a triangle normal gives a vector
on the triangle which is perpendicular to the edge, which is at an angle
$\frac{\Phi}{2}$ from the edge normal.

The full vertex shape tensor is a sum of the edge tensor, weighted by the
match of the normal of the edge to the normal of the vertex.

$S\left(i\right)=\frac{1}{A\left(i\right)}\sum_{j}\hat{N}\left(i\right)\cdot\hat{N}_{ij}\;S_{ij}$
(17)

We then project this 3x3 matrix in real space $\hat{x},\hat{y},\hat{z}$ to the
tangent plane of the vertex $\hat{d},\hat{t}$

$S_{2\times 2}=\begin{pmatrix}\hat{d}\cdot S\cdot\hat{d}&\hat{d}\cdot
S\cdot\hat{t}\\\ \hat{t}\cdot S\cdot\hat{d}&\hat{t}\cdot
S\cdot\hat{t}\end{pmatrix}$ (18)

The mean curvature $H$ at the vertex is half the trace, while the Gaussian
curvature $K$ is the determinant, which are two of the degrees of freedom in
the Hamiltonian.

The mismatch matrix can be calculated

$M=\begin{pmatrix}S_{dd}-\frac{C_{0}+D_{0}}{2}&S_{dt}\\\
S_{td}&S_{tt}-\frac{C_{0}-D_{0}}{2}\end{pmatrix}$ (19)

The angle between the director and the eigenvectors of the shape matrix is
what results in the $\omega$ angle dependence, which is the final degree of
freedom. The mismatch tensor is simply plugged to the equation 13 to compute
the bending energy of the vertex. If we integrate it across the whole
membrane, we get equation 1.

####  6.2.3 Dicretization of the order parameter

The order parameter is calculated by a python script shown below. The director
gives the direction of the inclusion $\vec{n}$.

⬇

import numpy as np

import sys

sys.path.append(’path_to_vtu_file’)

import vtu

v=vtu.PyVtu(’timestep_x.vtu’)

nvtx=len(v.type)

nematic_order=np.zeros(nvtx)

for i in np.arange(nvtx)[v.type!=4]:

js = v.get_neighbors(i)

js = js[v.c0[js]>0]

di = v.director[i]

nematic_order[i] = np.nanmean(3*(v.director[js]@di)**2-1)/2

v.update_array(’nematic_order’,nematic_order)

v.update_all()

v.write_vtu(’nematic_order.vtu’)

## References

  * [1] Jacob N Israelachvili.  Surface forces.  In The Handbook of Surface Imaging and Visualization, pages 793–816. CRC Press, 2022.
  * [2] V Kralj-Iglič, M Remškar, G Vidmar, M Fošnarič, and A Iglič.  Deviatoric elasticity as a possible physical mechanism explaining collapse of inorganic micro and nanotubes.  Physics Letters A, 296(2):151–155, 2002.
  * [3] Veronika Kralj-Iglič, Sasa Svetina, and B Žekž.  Shapes of bilayer vesicles with membrane embedded molecules.  European biophysics journal, 24:311–321, 1996.
  * [4] Veronika Kralj-Iglic.  Stability of membranous nanostructures: a possible key mechanism in cancer progression.  International journal of nanomedicine, pages 3579–3596, 2012.
  * [5] JB Fournier.  Nontopological saddle-splay and curvature instabilities from anisotropic membrane inclusions.  Physical review letters, 76(23):4436, 1996.
  * [6] VS Markin.  Lateral organization of membranes and cell shapes.  Biophysical Journal, 36(1):1–19, 1981.
  * [7] S Leibler.  Curvature instability in membranes.  Journal de Physique, 47(3):507–516, 1986.
  * [8] Veronika Kralj-Iglič, V Heinrich, Saša Svetina, and Boštjan Žekš.  Free energy of closed membrane with anisotropic inclusions.  The European Physical Journal B-Condensed Matter and Complex Systems, 10(1):5–8, 1999.
  * [9] Aleš Iglič, Blaž Babnik, Klemen Bohinc, Miha Fošnarič, Henry Hägerstrand, and Veronika Kralj-Iglič.  On the role of anisotropy of membrane constituents in formation of a membrane neck during budding of a multicomponent membrane.  Journal of Biomechanics, 40(3):579–585, 2007.
  * [10] Miha Fošnarič, Aleš Iglič, and Sylvio May.  Influence of rigid inclusions on the bending elasticity of a lipid membrane.  Physical Review E, 74(5):051503, 2006.
  * [11] Miha Fošnarič, Samo Penič, Aleš Iglič, Veronika Kralj-Iglič, Mitja Drab, and Nir Gov.  Theoretical study of vesicle shapes driven by coupling curved proteins and active cytoskeletal forces.  Soft Matter, 2019.
  * [12] NS Gov.  Guided by curvature: shaping cells by coupling curved membrane proteins and cytoskeletal forces.  Philosophical Transactions of the Royal Society B: Biological Sciences, 373(1747):20170115, 2018.
  * [13] Henry Hägerstrand, Lucyna Mrowczynska, Ulrich Salzer, Rainer Prohaska, Kimmo A Michelsen, Veronika Kralj-Iglič, and Aleš Iglič.  Curvature-dependent lateral distribution of raft markers in the human erythrocyte membrane.  Molecular Membrane Biology, 23(3):277–288, 2006.
  * [14] Aleš Iglič.  A possible mechanism determining the stability of spiculated red blood cells.  Journal of Biomechanics, 30(1):35–40, 1997.
  * [15] Arijit Mahapatra, Can Uysalel, and Padmini Rangamani.  The mechanics and thermodynamics of tubule formation in biological membranes.  The Journal of Membrane Biology, pages 1–19, 2021.
  * [16] Haleh Alimohamadi, Miriam K Bell, Shelley Halpain, and Padmini Rangamani.  Mechanical principles governing the shapes of dendritic spines.  Frontiers in Physiology, 12:657074, 2021.
  * [17] Aleš Iglič, Maruša Lokar, Blaž Babnik, Tomaž Slivnik, Peter Veranič, Henry Hägerstrand, and Veronika Kralj-Iglič.  Possible role of flexible red blood cell membrane nanodomains in the growth and stability of membrane nanotubes.  Blood Cells, Molecules, and Diseases, 39(1):14–23, 2007.
  * [18] Alex Veksler and Nir S Gov.  Phase transitions of the coupled membrane-cytoskeleton modify cellular shape.  Biophysical Journal, 93(11):3798–3810, 2007.
  * [19] Bojan Božič, Veronika Kralj-Iglič, and Saša Svetina.  Coupling between vesicle shape and lateral distribution of mobile membrane inclusions.  Physical Review E, 73(4):041915, 2006.
  * [20] Mitja Drab, David Stopar, Veronika Kralj-Iglič, and Aleš Iglič.  Inception mechanisms of tunneling nanotubes.  Cells, 8(6):626, 2019.
  * [21] Barbara Schamberger, Ricardo Ziege, Karine Anselme, Martine Ben Amar, Michał Bykowski, André PG Castro, Amaia Cipitria, Rhoslyn A Coles, Rumiana Dimova, Michaela Eder, et al.  Curvature in biological systems: Its quantification, emergence, and implications across the scales.  Advanced Materials, 35(13):2206110, 2023.
  * [22] Harvey T McMahon and Jennifer L Gallop.  Membrane curvature and mechanisms of dynamic cell membrane remodelling.  Nature, 438(7068):590–596, 2005.
  * [23] Joshua Zimmerberg and Michael M Kozlov.  How proteins produce cellular membrane curvature.  Nature reviews Molecular cell biology, 7(1):9–19, 2006.
  * [24] Mitja Drab, Žiga Pandur, Samo Penič, Aleš Iglič, Veronika Kralj-Iglič, and David Stopar.  A monte-carlo study of giant vesicle morphologies in nonequilibrium environments.  Biophysical Journal, 2021.
  * [25] Andrea Disanza, Sara Bisi, Moritz Winterhoff, Francesca Milanesi, Dmitry S Ushakov, David Kast, Paola Marighetti, Guillaume Romet-Lemonne, Hans-Michael Müller, Walter Nickel, et al.  Cdc42 switches irsp53 from inhibition of actin growth to elongation by clustering of vasp.  The EMBO journal, 32(20):2735–2750, 2013.
  * [26] Žiga Pandur, Samo Penič, Aleš Iglič, Veronika Kralj-Iglič, David Stopar, and Mitja Drab.  Surfactin molecules with a cone-like structure promote the formation of membrane domains with negative spontaneous curvature and induce membrane invaginations.  Journal of colloid and interface science, 650:1193–1200, 2023.
  * [27] Mitja Drab, Raj Kumar Sadhu, Yoav Ravid, Aleš Iglič, Veronika Kralj-Iglič, and Nir S Gov.  Modeling cellular shape changes in the presence of curved membrane proteins and active cytoskeletal forces.  Plasma Membrane Shaping, pages 415–429, 2023.
  * [28] Sho C Takatori and Amaresh Sahu.  Active contact forces drive nonequilibrium fluctuations in membrane vesicles.  Physical Review Letters, 124(15):158102, 2020.
  * [29] Gloria Mancinelli, Lucas Lamparter, Georgii Nosov, Tanumoy Saha, Anna Pawluchin, Rainer Kurre, Christiane Rasch, Mirsana Ebrahimkutty, Jürgen Klingauf, and Milos Galic.  Dendrite tapering actuates a self-organizing signaling circuit for stochastic filopodia initiation in neurons.  Proceedings of the National Academy of Sciences, 118(43):e2106921118, 2021.
  * [30] Francesco Cagnetta, Viktor Škultéty, Martin R Evans, and Davide Marenduzzo.  Universal properties of active membranes.  Physical Review E, 105(1):L012604, 2022.
  * [31] Ahis Shrestha, Fabien Pinaud, and Christoph A Haselwandter.  Mechanics of cup-shaped caveolae.  Physical Review E, 104(2):L022401, 2021.
  * [32] Ajoy Maji, Kinjal Dasbiswas, and Yitzhak Rabib.  Shape transitions in network model of active elastic shells.  arXiv preprint arXiv:2306.00735, 2023.
  * [33] Jeanne Stachowiak, Feng Yuan, Haleh Alimohamadi, Brandon Bakka, Andrea Trementozzi, Kasey Day, Nicolas Fawzi, and Padmini Rangamani.  Membrane bending by protein phase separation.  Proceedings of the National Academy of Sciences of USA, 118(11), 2021.
  * [34] Aleš Iglič, Blaž Babnik, Ulrike Gimsa, and Veronika Kralj-Iglič.  On the role of membrane anisotropy in the beading transition of undulated tubular membrane structures.  Journal of Physics A: Mathematical and General, 38(40):8527, 2005.
  * [35] Veronika Kralj-Iglič, V Heinrich, Saša Svetina, and Boštjan Žekš.  Free energy of closed membrane with anisotropic inclusions.  The European Physical Journal B-Condensed Matter and Complex Systems, 10:5–8, 1999.
  * [36] Wolfgang Helfrich.  Elastic properties of lipid bilayers: theory and possible experiments.  Zeitschrift für Naturforschung C, 28(11-12):693–703, 1973.
  * [37] Magali Deleu, Olivier Bouffioux, Hary Razafindralambo, Michel Paquot, Choukri Hbid, Philippe Thonart, Philippe Jacques, and Robert Brasseur.  Interaction of surfactin with membranes: a computational approach.  Langmuir, 19(8):3377–3385, 2003.
  * [38] Frederick C Frank.  I. liquid crystals. on the theory of liquid crystals.  Discussions of the Faraday Society, 25:19–28, 1958.
  * [39] Gordon Lasher.  Monte carlo results for a discrete-lattice model of nematic ordering.  Physical Review A, 5(3):1350, 1972.
  * [40] Netasan Ramakrishnan, PB Sunil Kumar, and John H Ipsen.  Monte carlo simulations of fluid vesicles with in-plane orientational ordering.  Physical Review E, 81(4):041922, 2010.
  * [41] Shuliang Chen, Tanvi Desai, James A McNew, Patrick Gerard, Peter J Novick, and Susan Ferro-Novick.  Lunapark stabilizes nascent three-way junctions in the endoplasmic reticulum.  Proceedings of the National Academy of Sciences, 112(2):418–423, 2015.
  * [42] Shuliang Chen, Peter Novick, and Susan Ferro-Novick.  Er network formation requires a balance of the dynamin-like gtpase sey1p and the lunapark family member lnp1p.  Nature cell biology, 14(7):707–716, 2012.
  * [43] Yoav Ravid, Samo Penič, Yuko Mimori-Kiyosue, Shiro Suetsugu, Aleš Iglič, and Nir S Gov.  Theoretical model of membrane protrusions driven by curved active proteins.  Frontiers in Molecular Biosciences, 10:1153420, 2023.
  * [44] G Gompper and DM Kroll.  Random surface discretizations and the renormalization of the bending rigidity.  Journal de Physique I, 6(10):1305–1320, 1996.
  * [45] G Gompper and DM Kroll.  Triangulated-surface models of fluctuating membranes.  Statistical Mechanics of Membranes and Surfaces, pages 359–426, 2004.
  * [46] Samo Penič, Aleš Iglič, Isak Bivas, and Miha Fošnarič.  Bending elasticity of vesicle membranes studied by monte carlo simulations of vesicle thermal shape fluctuations.  Soft Matter, 11:5004–5009, 2015.
  * [47] Miha Fošnarič, Samo Penič, Aleš Iglič, Veronika Kralj-Iglič, Mitja Drab, and Nir S Gov.  Theoretical study of vesicle shapes driven by coupling curved proteins and active cytoskeletal forces.  Soft Matter, 15(26):5319–5330, 2019.