Modelling fluid deformable surfaces with an emphasis on biological interfaces

TL;DR

This paper introduces a comprehensive continuum mechanics framework and computational methods for modeling fluid deformable surfaces, such as lipid bilayers and cell cortices, capturing their complex biological interactions and shape dynamics.

Contribution

It extends ALE formulations to deforming surfaces and develops a nonlinear Onsager formalism for coupled field and geometry modeling.

Findings

First 3D simulations of lipid bilayers and cell cortex models.

New computational methods for stiff PDE systems.

Framework captures finite curvature and shape changes.

Abstract

Fluid deformable surfaces are ubiquitous in cell and tissue biology, including lipid bilayers, the actomyosin cortex, or epithelial cell sheets. These interfaces exhibit a complex interplay between elasticity, low Reynolds number interfacial hydrodynamics, chemistry, and geometry, and govern important biological processes such as cellular traffic, division, migration, or tissue morphogenesis. To address the modelling challenges posed by this class of problems, in which interfacial phenomena tightly interact with the shape and dynamics of the surface, we develop a general continuum mechanics and computational framework for fluid deformable surfaces. The dual solid-fluid nature of fluid deformable surfaces challenges classical Lagrangian or Eulerian descriptions of deforming bodies. Here, we extend the notion of Arbitrarily Lagrangian-Eulerian (ALE) formulations, well-established for bulk media, to deforming surfaces. To systematically develop models for fluid deformable surfaces, which consistently treat all couplings between fields and geometry, we follow a nonlinear Onsager formalism according to which the dynamics minimize a Rayleighian functional where dissipation, power input and energy release rate compete. Finally, we propose new computational methods, which build on Onsager's formalism and our ALE formulation, to deal with the resulting stiff system of higher-order of partial differential equations. We apply our theoretical and computational methodology to classical models for lipid bilayers and the cell cortex. The methods developed here allow us to formulate/simulate these models for the first time in their full three-dimensional generality, accounting for finite curvatures and finite shape changes.

Equations

$$\bar{\mathbf{V}}\left(\mathbf{\xi },t\right)=\frac{d}{dt}{\mathbf{\varphi }}_{\mathbf{\xi }}(t)={\partial }_t\mathbf{\varphi }\left(\mathbf{\xi },t\right).$$

$$\mathbf{V}(\mathbf{x},t)=\bar{\mathbf{V}}\circ {\mathbf{\varphi }}_t^{-1}(\mathbf{x}),$$

$$\mathbf{V}=\mathbf{v}+v_n\mathbf{N},$$

$$\tilde{\mathbf{W}}(\mathbf{\xi },t)$$

$$\mathbf{W}(\mathbf{x},t)$$

$$\bar{\mathbf{c}}(\mathbf{\xi },t)$$

$$\tilde{\mathbf{c}}(\mathbf{\xi },t)$$

$$\mathbf{c}={\mathbf{\psi }}_{t\ast }\tilde{\mathbf{c}}=\left[D{\mathbf{\psi }}_t\tilde{\mathbf{c}}\right]\circ {\mathbf{\psi }}_t^{-1},$$

$$c^{\alpha }=\left({\partial }_b{\psi }_t^{\alpha }{\tilde{c}}^b\right)\circ {\mathbf{\psi }}_t^{-1},$$

$$c^a={\tilde{c}}^a\circ {\mathbf{\psi }}_t^{-1}.$$

$$\mathbf{V}=\mathbf{W}+\mathbf{c},$$

$$v_n=w_n,$$

$$\mathbf{v}=\mathbf{w}+\mathbf{c}.$$

$${\partial }_t\mathbf{\chi }\circ {\mathbf{\chi }}_t^{-1}=v_n\mathbf{N}.$$

$$\bar{\mathbf{t}}={\mathbf{\varphi }}_t^{\ast }\mathbf{t},\hat{\mathbf{t}}={\mathbf{\chi }}_t^{\ast }\mathbf{t},\tilde{\mathbf{t}}={\mathbf{\psi }}_t^{\ast }\mathbf{t},$$

$$D_{t}f(\mathbf{x},t)\equiv {\frac{d}{dt}f\left({\mathbf{\varphi }}_{\mathbf{\xi }}(t),t\right)|}_{\mathbf{\xi }={\mathbf{\varphi }}_t^{-1}(\mathbf{x})}.$$

$$D_{t}f(\mathbf{x},t)$$

$$={\partial }_t\left(f\circ {\mathbf{\varphi }}_t\right)\circ {\mathbf{\varphi }}_t^{-1}(\mathbf{x})={\mathbf{\varphi }}_{t\ast }{\partial }_t\left({\mathbf{\varphi }}_t^{\ast }f\right)(\mathbf{x}),$$

$$D_{t}f=L_{\mathbf{V}}f={\mathbf{\varphi }}_{t\ast }{\partial }_t\left({\mathbf{\varphi }}_t^{\ast }f\right).$$

$${\tilde{\partial }}_{t}f\equiv L_{\mathbf{W}}f={\mathbf{\psi }}_{t\ast }{\partial }_t\left({\mathbf{\psi }}_t^{\ast }f\right),$$

$${\partial }_{t}f\equiv L_{v_n\mathbf{N}}\mathbf{t}={\mathbf{\chi }}_{t\ast }{\partial }_t\left({\mathbf{\chi }}_t^{\ast }f\right).$$

$$D_{t}f={\tilde{\partial }}_{t}f+\text{bnabla}f\text{bcdot}\mathbf{c},$$

$$\mathbf{v}\text{bcdot}\mathbf{w}=\mathbf{g}\left(\mathbf{v},\mathbf{w}\right)=g_{ab}{v}^a{w}^b,$$

$$g_{ab}={\mathbf{e}}_a\text{bcdot}{\mathbf{e}}_b.$$

$$\ell [\mathbf{\alpha }]={\int }_{-1}^1|{\partial }_{\lambda }\mathbf{\alpha }|d\lambda ,$$

$$\cos \theta ={\left[\frac{{\partial }_{\lambda }\mathbf{\alpha }\text{bcdot}{\partial }_{\lambda }\mathbf{\beta }}{|{\partial }_{\lambda }\mathbf{\alpha }||{\partial }_{\lambda }\mathbf{\beta }|}\right]}_{\lambda =0}.$$

$${\left({\partial }_{\lambda }\mathbf{\alpha }\text{bcdot}{\partial }_{\lambda }\mathbf{\beta }\right)}_{\lambda =0}={\left[\left(\mathbf{g}\circ \mathbf{\alpha }\right)\left({\partial }_{\lambda }\mathbf{\alpha },{\partial }_{\lambda }\mathbf{\beta }\right)\right]}_{\lambda =0}={\left[\left(\bar{\mathbf{g}}\circ \bar{\mathbf{\alpha }}\right)\left({\bar{\mathbf{\alpha }}}^{\prime },{\bar{\mathbf{\beta }}}^{\prime }\right)\right]}_{\lambda =0}.$$

$$\left|{\partial }_{\lambda }\mathbf{\alpha }\right|=\sqrt{\left(\bar{\mathbf{g}}\circ \bar{\mathbf{\alpha }}\right)\left({\bar{\mathbf{\alpha }}}^{\prime },{\bar{\mathbf{\alpha }}}^{\prime }\right)}.$$

$$\bar{\mathbf{d}}=\frac{1}{2}{\partial }_t\bar{\mathbf{g}},$$

$$\mathbf{d}=\frac{1}{2}{\mathbf{\varphi }}_{t\ast }{\partial }_t\bar{\mathbf{g}}=\frac{1}{2}{\mathbf{\varphi }}_{t\ast }{\partial }_t\left({\mathbf{\varphi }}^{\ast }\mathbf{g}\right)=\frac{1}{2}{L}_{\mathbf{V}}\mathbf{g},$$

$${\frac{d}{dt}\left({\partial }_{\lambda }\mathbf{\alpha }\text{bcdot}{\partial }_{\lambda }\mathbf{\beta }\right)|}_{\lambda =0}=2{\left[\left(\mathbf{d}\circ \mathbf{\alpha }\right)\left({\partial }_{\lambda }\mathbf{\alpha },{\partial }_{\lambda }\mathbf{\beta }\right)\right]}_{\lambda =0},$$

$$\left[\frac{d}{dt}\left|{\partial }_{\lambda }\mathbf{\alpha }\right|\right]=\frac{1}{\left|{\partial }_{\lambda }\mathbf{\alpha }\right|}\left(\mathbf{d}\circ \mathbf{\alpha }\right)\left({\partial }_{\lambda }\mathbf{\alpha },{\partial }_{\lambda }\mathbf{\alpha }\right).$$

Full text

1 Introduction

Fluid deformable surfaces are a common motif in cell and tissue biology. For instance, lipid bilayers are fluid thin sheets that define the boundary of cells and compartmentalize them. They are the base material for the plasma membrane, the endoplasmic reticulum, mitochondria, or the Golgi apparatus. From a mechanical viewpoint, lipid bilayers are remarkable soft materials exhibiting a solid-fluid duality: while they store elastic energy when stretched or bent, as solid shells (Lipowsky, 1991), they cannot store elastic energy under in-plane shear, a situation under which they flow as viscous two-dimensional fluids (Dimova et al., 2006). This solid-fluid duality is tightly intertwined with membrane geometry: shape changes induce lipid flows that bring material from one part of the membrane to another (Evans & Yeung, 1994), whereas flows in the presence of curvature generate out-of-plane forces, which further curve the membrane (Rahimi et al., 2013). The solid-fluid duality of membranes is essential for cell function; it is required during cell motility and migration (Arroyo et al., 2012; Lieber et al., 2015), membrane trafficking (Sprong et al., 2001; Rustom et al., 2004), or to enable the mechano-adaptation of cells to stress (Staykova et al., 2013; Kosmalska et al., 2015). Furthermore, the in-plane fluidity of the membrane allows membrane inclusions, such as proteins, to diffuse (Sens et al., 2008). On the other hand, lipid bilayers are chemically responsive. Chemo-mechanical couplings can trigger tubulation (Roux et al., 2002), phase separation (Bacia et al., 2005), budding and fission (Staneva et al., 2004; Zhou & Yan, 2005), or pearling (Khalifat et al., 2014).

Another important instance of biological fluid surface is the actomyosin cortex, a thin network of cross-linked actin filaments lying immediately beneath the plasma membrane of animal cells (Bray & White, 1988). Within this network, myosin motors exert active forces by consuming chemical energy in the form of adenosine triphosphate (ATP), that generate active tension (Salbreux et al., 2012). Furthermore, the cell cortex undergoes dynamic remodelling, or turnover, in less than one minute, as a result of the polymerization and depolymerization of actin and the binding and unbinding of cross-linkers (Howard, 2001). The cell cortex behaves as an elastic network at short time-scales and as a quasi-two-dimensional viscous fluid at longer time-scales due to turnover. The interplay between remodelling, elasticity, and active forces in this thin cortical layer plays a critical role in different cellular processes such as cytokinesis (Levayer & Lecuit, 2012), or migration (Bergert et al., 2015; Ruprecht et al., 2015; Callan-Jones et al., 2016), where the coupling between shape and actin flows becomes apparent.

In summary, fluid deformable surfaces are ubiquitous interfaces in biology, adopting three dimensional dynamical shapes, involving chemo-mechanical couplings, and exhibiting a dual solid-fluid behaviour. The mechanics of these biological interfaces plays an essential role in processes from the subcelullar to the tissue scale (not discussed here). However and despite recent efforts (Salbreux & Jülicher, 2017; Sahu et al., 2017; Sauer et al., 2017), a general theoretical and computational framework to describe the multiphysics and geometry-dependent mechanics of these systems has been lacking. Towards filling this gap, here we develop a three-dimensional non-linear modelling and simulation framework for fluid deformable surfaces. Even though such interfacial fluids are embedded in a bulk fluid or confined to substrates, here we focus only on the mechanics of the surface. The coupling of interfaces with a bulk fluid (Salac & Miksis, 2011; Woodhouse & Goldstein, 2012; Farutin & Misbah, 2012; Shen et al., 2018; Laadhari et al., 2017) or a substrate (Staykova et al., 2013) has been examined extensively in other works.

Different mathematical representations of the kinematics of fluid deformable surfaces have been proposed. In a common approach (Secomb & Skalak, 1982; Barthes-Biesel & Sgaier, 1985), the velocity field of the interface is defined as the restriction to the surface of the velocity field of the bulk fluid assuming a no-slip condition. The governing equations for the interface are then obtained in terms of time-dependent projection operators. Even if the interaction with a bulk fluid is not considered, the interfacial velocity can be extended to a tubular neighbourhood in 3D around the interface to find the governing equations, which can be shown to be independent of the extension (Dziuk & Elliott, 2007, 2013). This approach has been applied to the numerical simulations of lipid membranes (Rodrigues et al., 2015; Barrett et al., 2015, 2016a) and to the numerical solution of vector PDEs, such as Navier-Stokes, on surfaces of fixed shape (Hansbo et al., 2016; Reuther & Voigt, 2018; Fries, 2018). This framework can be adapted to parametrization-free descriptions of surfaces using level-sets (Dziuk & Elliott, 2013; Burman et al., 2015). However, by extending the problem to Euclidean space, this approach hides much of the geometric structure of the governing equations. Furthermore, it is not obvious how to extend it to bilayer interfaces such as lipid membranes, in which individual monolayers are bound to the mid-surface but can slip relative to each other. An alternative approach, pioneered by Scriven (1960), distinguishes between the intrinsic (tangential) velocity of particles as seen from within the surface, and the extrinsic (normal) surface velocity, which changes its shape and thus its metric tensor (Aris, 1962). This approach, revisited in different theoretical and computational works (Hu et al., 2007; Arroyo & Desimone, 2009; Rahimi & Arroyo, 2012; Sahu et al., 2017), requires the language and computational tools of differential geometry, provides a clear geometric picture of the governing equations, and eloquently shows the tight interplay between shape changes and interfacial flows. Here, we show that, by decoupling shape changes and tangential flows, this approach can naturally generalize Arbitrarily Lagrangian-Eulerian (ALE) methods, well established for bulk media (Hirt et al., 1974; Donea & Huerta, 2003), to fluid deformable surfaces. Thus, this formalism (1) alleviates the large distortions of a pure Lagrangian framework, which usually requires intensive remeshing (Rodrigues et al., 2015), and (2) allows us to deal with multilayer systems by considering independent tangential velocities for each monolayer.

To deal with the multiphysics aspects of fluid surfaces, we base our approach on a nonlinear Onsager’s formalism (Arroyo et al., 2018; Doi, 2011; Mielke, 2012; Peletier, 2014), which provides a unified variational framework for the dissipative dynamics of soft-matter systems. In this formalism, the dynamics minimize a Rayleighian functional and result from the interplay between energetic driving forces, dissipative drag forces and external forces, each of them deriving from potentials that are the sum of individual contributions for each physical mechanism. Complex models coupling different physics can be assembled by just adding more terms to the energy and dissipation potentials, and encoding in them the interactions between the different physical mechanisms. Thus, this framework provides a transparent and thermodynamically consistent method to generate complex models. Onsager’s formalism is applicable to capillarity, elasticity, low Reynolds number hydrodynamics, reaction-diffusion systems, and provides a natural framework to model biological activity. In different contexts, similar ideas have been referred to by different names, such as extremal principles in non-equilibrium thermodynamics studied in physics (Martyushev & Seleznev, 2006; Lebon et al., 2008), in materials modelling (Ziegler, 1958; Ziegler & Wehrli, 1987; Ortiz & Stainier, 1999; Fischer et al., 2014) or in atmospheric transport processes (Paltridge, 1975). The Onsager formalism used here generalizes previous minimum principles identified in low Reynolds number hydrodynamics coupled to capillary (Skalak, 1970) or viscoelastic interfaces (Secomb & Skalak, 1982; Dörries & Foltin, 1996).

In addition to the geometric and multiphysics aspects of the theory, the three-dimensional simulation of fluid surfaces requires specialized numerical methods since the resulting equations (1) usually involve higher-order derivatives of the parametrization, (2) lead to a mixed system of elliptic and hyperbolic partial differential equations and (3) are stiff and difficult to integrate in time (Rahimi et al., 2013). Indeed, surface shape enters into the energy and dissipation expressions through curvature, which involves second-order derivatives of the parametrization. From a finite element method (FEM) perspective, this requires the basis functions parametrizing the surface to be in $H^{2}$ (square-integrable functions whose first- and second-order derivatives are also square-integrable). Here, we resort to subdivision surfaces, which have already been used to study the equilibrium shapes of lipid bilayers (Feng & Klug, 2006; Ma & Klug, 2008) and to analyze thin shells (Cirak et al., 2000; Cirak & Ortiz, 2001; Cirak & Long, 2011; Zhang & Arroyo, 2014; Li et al., 2018). Based on a time-incremental version of Onsager’s formalism, we develop variational time-integrators (Ortiz & Stainier, 1999; Peco et al., 2013), which are nonlinearly and unconditionally stable and allow us to adapt the time-step spanning orders of magnitude during the dynamics of fluid deformable surfaces.

The paper is structured as follows. In section 2, we develop a theoretical description of fluid surfaces, including Lagrangian, Eulerian and ALE formulations. We introduce the rate-of-deformation tensor and the Reynolds transport theorem. We also describe a useful set of tools to represent the kinematics of fluid deformable surfaces. In section 3, we describe several classical models of fluid surfaces to describe the dynamics of lipid bilayers and the cell cortex. We show how Onsager’s variational formalism provides a direct and transparent tool to derive complex governing equations. In section 4, we describe the discretization, both in time and space, of the equations governing the dynamics of general fluid surfaces. We introduce a variational time-integrator based on Onsager’s formalism, and show how to discretize the different fields defined on the surface. In section 5, we exercise the models in section 3 through several examples simulated using the techniques described in section 4. Finally, we conclude in section 6 with a summary and discussion of the manuscript, along with suggestions for future work.

2 Mathematical description of fluid deformable surfaces

In this section, we mathematically describe fluid surfaces as a two-dimensional continua moving and deforming in Euclidean space. One way to represent this kind of systems is through a Lagrangian parametrization of the surface, $\bm{\phi}(\bm{\xi},t)$, in which a material particle is identified with a point $\bm{\xi}^{*}$ in parametric domain and $\bm{\phi}(\bm{\xi}^{*},t)$ follows its trajectory in time. However, Lagrangian parametrizations present two major drawbacks for the description of fluid surfaces. First, due to the fluid nature of the interface, Lagrangian parametrizations suffer from very large distortions that are difficult to accommodate with conventional discretization schemes. Second, a single Lagrangian parametrization cannot track simultaneously all material particles in a multilayer interface. For example, in a lipid bilayer, two material particles representing lipid molecules from each monolayer occupy the same position on the surface. A single Lagrangian parametrization cannot track the time-evolution of both simultaneously because they can slip relative to each other.

In this section, we examine the definition of Lagrangian, Eulerian and ALE parametrizations of material surfaces and establish their relations. Associated with the flow generated by these parametrizations, we define the Lagrangian, Eulerian and ALE time-derivatives of fields on the surface. We then introduce the right Cauchy deformation tensor and the rate-of-deformation tensor, which characterizes the rate at which lengths, angles and areas transform on the time-evolving surface. We examine time-derivatives of integrals on time-evolving surfaces, and derive the form of Reynolds transport theorem and conservation of mass for the Lagrangian, Eulerian and ALE descriptions. Finally, we introduce some mathematical tools to represent the kinematics of fluid surfaces.

Throughout the manuscript, we make extensive use of the differential geometry of surfaces, including the definition of the metric tensor or first fundamental form $\bm{g}$, second fundamental form or shape operator $\bm{k}$, covariant differentiation $\bnabla$, and Lie derivation $L_{\bm{v}}$, along with push-forwards and pull-backs by maps. Contravariant components of a tensor are denoted by superscripts, whereas covariant components are denoted by subscripts; for instance, the components of the metric tensor are denoted by $g_{ab}$, whereas the components of a tangent vector are denoted by $v^{a}$. We use Latin letters to denote indices running from $1$ to $2$, representing tensors on the tangent space of the surface, and Greek letters to denote indices running from $1$ to $3$, used for tensors in Euclidean space. We follow Einstein’s notation: contravariant and covariant indices with the same label are implicitly summed $T_{a\cdots}T^{a\cdots}=\sum_{a=1}^{2}T_{a\cdots}T^{a\cdots}$. We refer to do Carmo (2016); Do Carmo (1992); Willmore (1996) for background texts about the differential geometry of surfaces and manifolds.

2.1 Lagrangian, Eulerian and ALE parametrizations

We consider the parametrization of a two-dimensional continuum $\Gamma_{t}$ moving and deforming in $\mathbb{R}^{3}$. In a Lagrangian parametrization of $\Gamma_{t}$, $\bm{\phi}:\bar{\Gamma}\times\mathcal{I}\ni\left(\bm{\xi},t\right)\mapsto\bm{x}\in\Gamma_{t}$, where $\bar{\Gamma}\subset\mathbb{R}^{2}$ and $\mathcal{I}\subset\mathbb{R}$, a point $\bm{\xi}=(\xi_{1},\xi_{2})\in\bar{\Gamma}$ identifies a material particle and the curve obtained by fixing $\bm{\xi}$, $\bm{\phi}_{\bm{\xi}}(t)=\bm{\phi}\left(\bm{\xi},t\right)$, is its trajectory in $\mathbb{R}^{3}$ (see figure 1). We focus on a specific chart, although the arguments presented in this section can be trivially extended to surfaces covered by an atlas of charts. For systems with multiple components, e.g. multilayer systems, where material particles of different components coexist at the same point $\bm{x}\in\Gamma_{t}$, a single Lagrangian parametrization of $\Gamma_{t}$ does not exist. This is the case of a lipid bilayer, where a point $\bm{x}\in\Gamma_{t}$ has simultaneously attached two material particles belonging to each monolayer. Nevertheless, we can always define a Lagrangian parametrization for each of the components of the system independently so that the results in this and following sections can be applied to each component (monolayer) separately. The time derivative of the Lagrangian parametrization is the material velocity

[TABLE]

The spatial velocity $\bm{V}$ on $\Gamma_{t}$ is obtained by composition with $\bm{\phi}^{-1}_{t}$

[TABLE]

where $\bm{\phi}_{t}(\bm{\xi})=\bm{\phi}(\bm{\xi},t)$ is obtained by fixing time $t$. In general, $\bm{V}$ has a tangential and a normal component to $\Gamma_{t}$

[TABLE]

where $\bm{N}$ is the unit normal to the surface. The normal velocity $v_{n}$ characterizes shape changes of $\Gamma_{t}$ while $\bm{v}$ represents the flow of material tangent to $\Gamma_{t}$. In the remainder of the paper we denote by upper-case letters vectors with tangential and normal components to $\Gamma_{t}$ and by lower-case letters vectors that are tangent to $\Gamma_{t}$. We now introduce an alternative parametrization of the surface $\bm{\psi}:\tilde{\Gamma}\times\mathcal{I}\ni(\bm{\xi},t)\mapsto\bm{x}\in\Gamma_{t}$, where $\tilde{\Gamma}\subset\mathbb{R}^{2}$. The curves of constant $\bm{\xi}$, $\bm{\psi}_{\bm{\xi}}(t)=\bm{\psi}(\bm{\xi},t)$ do not follow trajectories of material particles in general.

The velocity fields associated with this parametrization are

[TABLE]

We can construct a map relating both parametrizations $\bm{\theta}=\bm{\psi}^{-1}_{t}\circ\bm{\phi}:\bar{\Gamma}\times\mathcal{I}\rightarrow\tilde{\Gamma}$, a diffeomorphism from $\mathbb{R}^{2}$ to $\mathbb{R}^{2}$ at each $t$. The curves of constant $\bm{\xi}$, $\bm{\theta}_{\bm{\xi}}(t)=\bm{\theta}(\bm{\xi},t)$, track the parametric positions of material particles evolving in $\tilde{\Gamma}$, and have a velocity

[TABLE]

To physically interpret $\tilde{\bm{c}}$, we define its push-forward by $\bm{\psi}_{t}$ as

[TABLE]

where $\bm{\psi}_{t*}$ denotes the push-forward, and $D\bm{\psi}_{t}$ stands for the differential of $\bm{\psi}_{t}$, a linear mapping from the tangent space of $\tilde{\Gamma}$ at $\bm{\xi}$, $T_{\bm{\xi}}\tilde{\Gamma}$, to the tangent space of $\Gamma_{t}$ at $\bm{x}=\bm{\psi}_{t}(\bm{\xi})$, $T_{\bm{x}}\Gamma_{t}$. The components of this vector in the global basis of Euclidean space are

[TABLE]

where we have used the notation $\partial_{a}=\partial_{\xi_{a}}$. This expression shows that that in the basis $\bm{e}_{a}=\partial_{a}\bm{\psi}_{t}\circ\bm{\psi}^{-1}_{t}(\bm{x})$ of the tangent space $T_{\bm{x}}\Gamma_{t}$, the convected basis by $\bm{\psi}$, the components of $\bm{c}$ are simply

[TABLE]

Thus, in the convected basis, the components in $T_{\bm{\psi}_{t}^{-1}(\bm{x)}}\tilde{\Gamma}$ coincide with those in $T_{\bm{x}}\Gamma_{t}$. Using the chain rule and previous definitions (see appendix A), we recover the classical relation between Lagrangian and ALE parametrizations in the bulk (Donea & Huerta, 2003),

[TABLE]

and thus $\bm{c}$ represents the relative velocity of material particles with respect to the parametrization given by $\bm{\psi}$. Since $\bm{c}$ is the push-forward of a vector field with respect to $\bm{\psi}$, then it is tangent to $\Gamma_{t}$. Comparing Eqs. (3) and (4), we conclude that

[TABLE]

and

[TABLE]

This reflects that, since both parametrizations describe the same shape, their normal velocities, characterizing shape changes, must coincide. With this in mind, we can now introduce the notion of Eulerian parametrization in the context of a time-evolving surface. We say that a parametrization $\bm{\chi}$ is Eulerian if its velocity field is always perpendicular to the surface

[TABLE]

In summary, the parametrization $\bm{\phi}$ is a Lagrangian parametrization that tracks the evolution of material particles as they move with and along $\Gamma_{t}$. On the other hand, $\bm{\chi}$ is an Eulerian parametrization whose velocity is always perpendicular to $\Gamma_{t}$ regardless of the tangential flows of material. These parametrizations are special cases of a general parametrization $\bm{\psi}$, which may present tangential movements not consistent with the velocity of material particles. This kind of parametrization is referred to as an arbitrarily Lagrangian-Eulerian (ALE) parametrization.

We introduce here some notation. We denote the pull-backs of a tensor $\bm{t}$ on $\Gamma_{t}$ by the Lagrangian, Eulerian and ALE maps by

[TABLE]

where $\bm{\phi}^{*}_{t}$ denotes the pull-back through $\bm{\phi}_{t}$.

2.2 Material, Eulerian and ALE time derivatives

We introduce next the concept of time-derivative of fields on $\Gamma_{t}$. Let us focus for simplicity on a scalar field over $\Gamma_{t}$, $f(\bm{x},t)$. We first note that the operator $\partial_{t}$ acting on $f(\bm{x},t)$, with the usual meaning of taking the time-derivative at fixed $\bm{x}$, is not well defined since $\bm{x}$ cannot be held fixed on a time-evolving surface in general (Cermelli et al., 2005). The idea of time-derivative can be more easily rationalized resorting to a parametrization. Let us first consider the Lagrangian parametrization $\bm{\phi}$. Fixing a point $\bm{\xi}\in\bar{\Gamma}$, we can compute how $f(\bm{x},t)$ changes along the curve $\bm{\phi}_{\bm{\xi}}(t)$. We define the material time derivative of the scalar $f$ as

[TABLE]

We note that $f\left(\bm{\phi}_{\bm{\xi}}(t),t\right)$ is a function of $t$ only and therefore the right-hand side of the previous expression is the usual derivative of a function of one variable. By noting that the pull-back of $f$ onto $\bar{\Gamma}$ is $\bar{f}=\bm{\phi}^{*}_{t}f=f\circ\bm{\phi}_{t}$, we can rewrite the previous expression as

[TABLE]

where $\partial_{t}\bar{f}$ has the usual meaning of taking the partial derivative of $\bar{f}$ at fixed $\bm{\xi}$. The last expression in this equation can be worded as the push-forward of the time-derivative of the pull-back of $f$ by the Lagrangian parametrization $\bm{\phi}$. This is nothing but the Lie-derivative of $f$ along the flow generated by $\bm{V}$, usually denoted by $L_{\bm{V}}f$, which is an extension to non-tangential vector fields of the usual definition of Lie-derivative (Do Carmo, 1992; Arroyo & Desimone, 2009). We note that the Lie-derivative depends on $\bm{\phi}$ only through $\bm{V}$. Thus, we can write the material time-derivative as

[TABLE]

We can equivalently define the ALE time-derivative of $f$

[TABLE]

and the Eulerian time derivative

[TABLE]

In the left-hand-side of this equation. the meaning of the symbol $\partial_{t}$ is clear: it measures the rate of change of $f$ along the flow normal to $\Gamma_{t}$. If the shape of $\Gamma_{t}$ remains stationary, then $\partial_{t}$ recovers the usual meaning of taking the derivative with respect to time at fixed $\bm{x}$. We note that the symbol $\partial_{t}$ retains the usual meaning when applied to fields on parametric domains, e.g. $\partial_{t}\bar{f}=\lim_{\Delta t\rightarrow 0}\left(\bar{f}(\bm{\xi},t+\Delta t)-\bar{f}(\bm{\xi},t)\right)/\Delta t$, and should not be confused with the definition equation (18) for fields on $\Gamma_{t}$. The operators $D_{t}$, $\tilde{\partial}_{t}$ and $\partial_{t}$ are related. For instance, using previous definitions in equation (13) and the chain rule (see appendix B), we find that

[TABLE]

where $\bnabla$ denotes the covariant derivative, or here the surface gradient, of $f$.

2.3 Rate-of-deformation tensor

An important tensor on $\Gamma_{t}$ is the first fundamental form or metric tensor $\bm{g}$. The metric tensor induces a scalar product on the tangent space of $\Gamma_{t}$ that allows us to measure lengths, angles and areas on $\Gamma_{t}$. Given two tangent vectors to $\Gamma_{t}$, $\bm{v}$ and $\bm{w}$, the scalar product is defined by

[TABLE]

where the notation $\bm{g}(\bcdot,\bcdot)$ views $\bm{g}$ as a bilinear form. For surfaces in $\mathbb{R}^{3}$, the metric tensor is defined so that the scalar product on $\Gamma_{t}$ coincides with the scalar product in $\mathbb{R}^{3}$. Then, given a basis $\{\bm{e}_{1},\bm{e}_{2}\}$ of $T\Gamma_{t}$, the tangent bundle of $\Gamma_{t}$, the components of the metric tensor in this basis are

[TABLE]

Let us consider two curves in the parametric domain $\bar{\Gamma}$, given by $\bar{\bm{\alpha}}(\lambda):[-1,1]\rightarrow\bar{\Gamma}$ and $\bar{\bm{\beta}}(\lambda):[-1,1]\rightarrow\bar{\Gamma}$, that cross at $\lambda=0$, and the image of these curves by $\bm{\phi}$, $\bm{\alpha}(\lambda,t)=\bm{\phi}(\bar{\bm{\alpha}}(\lambda),t)$ and $\bm{\beta}(\lambda,t)=\bm{\phi}(\bar{\bm{\beta}}(\lambda),t)$ (see figure 2).

The length of $\bm{\alpha}$ (and equivalently of $\bm{\beta}$) is given by the functional

[TABLE]

where $|\bm{v}|=\sqrt{\bm{v}\bcdot\bm{v}}$ is the norm of $\bm{v}$. The angle $\theta$ between curves $\bm{\alpha}$ and $\bm{\beta}$ at their point of intersection is given by

[TABLE]

The time-evolution of the lengths of material curves and the angles between them measures how the material deforms. It is interesting to note that the pull-back of $\bm{g}$, $\bar{\bm{g}}=\bm{\phi}^{*}\bm{g}$, induces a time-dependent scalar product on $\bar{\Gamma}$ that allows us to compute products of deformed vectors from their time-independent description on $\bar{\Gamma}$. For instance, one can easily see that

[TABLE]

Equivalently,

[TABLE]

Thus, scalar products, lengths and angles of material curves on $\Gamma_{t}$, such as $\bm{\alpha}$ and $\bm{\beta}$, can be measured on $\bar{\Gamma}$, from the time-independent $\bar{\bm{\alpha}}$ and $\bar{\bm{\beta}}$, with the time-dependent scalar product induced by $\bar{\bm{g}}$. It is clear from Eqs. (24) and (25) that the time-dependence of these measures of local deformation is completely encoded in $\bar{\bm{g}}$. We conclude that the tensor $\bar{\bm{g}}$ characterizes the deformation of $\Gamma_{t}$. In continuum mechanics, this tensor is referred to as the (right Cauchy-Green) deformation tensor and is generally denoted by $\bm{C}$. The time-derivative of this tensor defines a new tensor over $\bar{\Gamma}$

[TABLE]

where the $1/2$ is introduced here to follow the usual convention. The push-forward of this tensor to $\Gamma_{t}$ by $\bm{\phi}$ defines the so-called rate-of-deformation tensor,

[TABLE]

where we recognize again the structure of a Lie derivative, this time applied to the metric tensor. The rate of change of the scalar product can then be written as

[TABLE]

and the rate of change of the norm as

[TABLE]

Thus, the rate of change of local deformation of $\Gamma_{t}$ is encoded in $\bm{d}$. As shown in appendix C, see also Marsden & Hughes (1994); Z. Wu et al. (2005), the rate-of-deformation tensor for a surface moving in Euclidean space can be written as

[TABLE]

where $\bnabla$ is the covariant derivative, and $\bm{k}$ is the shape operator characterizing the local curvature of the surface and defined as

[TABLE]

From this expression, it is clear that the surface $\Gamma_{t}$ deforms through tangential flows, which contribute to the rate-of-deformation tensor with the usual term $\left[\bnabla\bm{v}+\left(\bnabla\bm{v}\right)^{T}\right]/2$, but also through the change in shape of $\Gamma_{t}$, which contributes with the term $-v_{n}\bm{k}$. This relation illustrates the coupling between tangential flows and shape changes in the presence of curvature.

2.4 Reynolds transport theorem and conservation of mass

In this section we extend the concept of Lagrangian, Eulerian and ALE time-derivatives of integrals on $\Gamma_{t}$. Consider a subset $\Xi\subset\Gamma_{t}$, a scalar field $f:\Gamma_{t}\rightarrow\mathbb{R}$, and define

[TABLE]

To compute this integral, we can pull-back $fdS$ onto $\bar{\Gamma}$

[TABLE]

where $\bar{\Xi}=\bm{\phi}_{t}^{-1}(\Xi)$, $\bar{J}=\sqrt{\bar{g}}$, $\bar{g}=\det(\bar{\bm{g}})=\det\left(D\bm{\phi}^{T}D\bm{\phi}\right)$ and $d\bm{\xi}=d\xi_{1}d\xi_{2}$. We define the material time derivative of $I$ as

[TABLE]

This expression characterizes the rate of change of the integral $I$ when the domain $\Xi$ is a material subset of $\Gamma_{t}$, i.e. it evolves following the flow generated by $\bm{\phi}$ (see figure 3).

Developing the definition, we have

[TABLE]

The rate of change of $\bar{J}=\sqrt{\bar{g}}$ can be written in terms of $\bm{d}$ by noting that $\partial_{t}\bar{J}=\partial_{t}\bar{g}/(2\bar{J})$ and using Jacobi’s formula $\partial_{t}\bar{g}=\bar{g}\bar{\bm{g}}^{-1}\bm{:}\left(\partial_{t}\bar{\bm{g}}\right)=\bar{g}\left[\bm{g}^{-1}\bm{:}\left(L_{\bm{V}}\bm{g}\right)\circ\bm{\phi}_{t}\right]=2\bar{g}\left[\text{tr}\bm{d}\circ\bm{\phi}_{t}\right]$, where $\text{tr}\bm{d}=d^{a}_{~{}a}=g^{ab}d_{ab}$ is the trace $\bm{d}$. Thus, we have

[TABLE]

where we have used equation (30), $\bnabla\bcdot\bm{v}=\bnabla_{a}v^{a}$ is the surface divergence of the tangential vector field $\bm{v}$, and we define the mean curvature as $H=g^{ab}k_{ab}$. Then,

[TABLE]

Using Eqs. (19) and the divergence theorem for surfaces, we can rewrite the previous equation in different ways

[TABLE]

where $\partial\Xi$ indicates the boundary curve of $\Xi$ and $\bm{m}$ the outer normal to $\partial\Xi$ and tangent to $\Xi$. Eqs. (37)-(41) are the equivalent to Reynold’s transport theorem for material domains in terms of the material, Eulerian and ALE time-derivative of $f$. As for scalar fields, we can extend the notion of material time-derivative of an integral relative to other parametrizations. In particular, we can consider the parametric domain $\tilde{\Xi}=\bm{\psi}_{t}^{-1}(\Xi)$, and the time-derivative

[TABLE]

where $\tilde{J}=\det\tilde{\bm{g}}$. This time-derivative characterizes the rate of change of $I$ when it follows the flow generated by the ALE parametrization. One can easily prove that

[TABLE]

For an Eulerian parametrization, one equivalently finds

[TABLE]

From these expressions, it is clear that for a closed surface $D_{t}I=\tilde{\partial}_{t}I=\partial_{t}I$.

The previous expression can be used to derive the statement of conservation of mass on fluid surfaces. Indeed, in the special case of $f=\rho$, the mass density per unit area, conservation of mass for every material sub-domain $\Gamma_{t}$ requires that

[TABLE]

where $r$ is the rate of creation of mass per unit area, which may for instance result from the exchange of material with the bulk. Since this must hold for every subdomain $\Xi$, we can localize the statement to obtain Lagrangian, Eulerian and ALE forms of local conservation of mass

[TABLE]

For inextensible fluid surfaces in the absence of mass exchange, balance of mass reduces to $D_{t}\rho=0$, leading to the condition

[TABLE]

Thus, for an inextensible surface with curvature, any shape change must be accompanied by a tangent flow to fulfill the inextensibility constraint, further highlighting the tight coupling between tangent flows and shape changes in the presence of curvature.

2.5 Representation of kinematics for fluid deformable surfaces

In previous sections, we have seen that Lagrangian parametrizations are natural tools to define the deformation tensor $\bar{\bm{g}}$, the rate-of-deformation tensor $\bm{d}$, and to establish the transport theorem on a time-deforming surface. A time-dependent Lagrangian parametrization contains information about shape changes ($v_{n}$) and about interfacial flows ($\bm{v}$). In practical computations, however, Lagrangian parametrizations are not well-suited for fluid surfaces because they exhibit large distortions, requiring intensive remeshing (Rodrigues et al., 2015), and because a single Lagrangian parametrization cannot describe a multicomponent system like a lipid bilayer, where monolayers can slip relative to each other. In this case, one could consider a Lagrangian parametrization for each component, which, however, increases the number of degrees of freedom since each parametrization describes both tangential motion and shape, whereas only tangential motions are independent of each other. In the present section, we provide a set of modelling tools, which are useful for a clean formulation of physical models of fluid surfaces and particularly for their numerical discretization.

In the previous section, we have introduced the notion of a time-dependent ALE parametrization $\bm{\psi}$ to describe the time-evolution of a material surface $\Gamma_{t}$, which can alleviate mesh distortion when dealing with fluid surfaces since it does not follow material particles. We note, however, that $\bm{\psi}$ does not contain information about the tangential motion of material particles (the interfacial flows) given by $\bm{v}$, since $\bm{v}$ and the tangential velocity of $\bm{\psi}$ differ by a relative velocity $\bm{c}$. This fact confronts us with two issues. First, how to select the tangential velocity of $\bm{\psi}$, which is arbitrary in the sense of not being prescribed by any physical law. Second, since $\bm{v}$ needs to be considered as an object independent of $\bm{\psi}$, how to parametrize tangential vector fields? The first issue has been addressed by introducing a numerical drag, which limits the tangential motion of $\bm{\psi}$ (Rahimi & Arroyo, 2012; Ma & Klug, 2008). One could also use the physically unconstrained tangential degrees of freedom of $\bm{\psi}$ to perform dynamical mesh adaptation (Veerapaneni et al., 2011). These approaches, however, increase the number of essential degrees of freedom required to describe shape changes (from one to three) and require parameter tuning. Instead, in section 2.5.1 we develop a special kind of ALE parametrization based on an offset (Rangarajan & Gao, 2015), which parametrizes $\bm{\psi}$ using a scalar field over $\Gamma_{t}$. Regarding the second issue, we note that interpolating tangent vector fields in a system with multiple charts is delicate, see section 4. In section 2.5.2, we introduce the Hodge decomposition of vector fields in terms of scalar fields, whose interpolation is straightforward.

2.5.1 An ALE parametrization based on an offset

We define next a restricted ALE parametrization, which by construction is devoid of the arbitrary freedom associated with tangential motions. Let us consider the state of the surface at a given time $t_{0}$, $\Gamma_{t_{0}}$, and a parametrization of this surface $\bm{\psi}_{0}(\bm{\xi})$. We consider a vector field $\bm{M}(\bm{\xi})$, representing a field of directors over $\Gamma_{t_{0}}$, with non-zero normal component but not necessarily coinciding with the normal field of $\Gamma_{t_{0}}$. We define a family of parametrizations of $\Gamma_{t}$ at time $t>t_{0}$ in terms of the offset of a point $\bm{x}=\bm{\psi}_{0}(\bm{\xi})$ along $\bm{M}(\bm{\xi})$,

[TABLE]

see figure 4A. The field that characterizes the time-evolution of the parametrization is $h$, a simple scalar field on $\tilde{\Gamma}$.

In principle, this approach is not completely general, since it only allows us to parametrize the surfaces lying in the so-called tubular neighbourhood of $\Gamma_{t_{0}}$ (do Carmo, 2016). For some interval $I=(t_{0},t_{0}+\delta t)$, the deformed surface $\Gamma_{t}$ will lie in the tubular neighbourhood of $\Gamma_{t_{0}}$ if the time-evolution is smooth. However, after some time, $\Gamma_{t}$ may leave the tubular neighbourhood of $\Gamma_{t_{0}}$ (see figure 4B for an example). A simple solution to this issue is to then update the reference configuration $\Gamma_{t_{0}}$. This kind of parametrization, proposed by Rangarajan & Gao (2015), generalizes the classical Monge parametrization, which is recovered by setting $\Gamma_{t_{0}}$ to a plane, $\bm{M}$ to its constant normal and $h$ to the height of the surface $\Gamma_{t}$ with respect to the plane (do Carmo, 2016). We finally note that for this kind of surface parametrization, we have

[TABLE]

Since $h$ is a scalar field on $\tilde{\Gamma}$, $\partial_{t}h$ in this equation has the usual meaning of time differentiation at fixed $\bm{\xi}$. In practice, $\bm{M}$ can be chosen to be $\bm{N}_{0}$, the field of normals of the reference surface as in (Rangarajan & Gao, 2015). This leads to an Eulerian parametrization at $t=t_{0}$, and close to it at later times. Thus, $\bm{W}$ will have in general non-zero tangential components, and therefore this parametrization is neither Eulerian nor Lagrangian. Instead, it is an ALE parametrization depending on a generalized height field $h$, in which the arbitrariness is removed by following equation (48) and choosing the field of directors $\bm{M}$.

2.5.2 Velocity potentials: Hodge decomposition

Given a vector field $\bm{V}\in\mathbb{R}^{3}$, it is well-known that $\bm{V}$ admits a decomposition in terms of the gradient of a function $\Phi$ and the curl of a vector potential $\bm{A}$ in what is called the Helmholtz decomposition,

[TABLE]

where here $\bnabla$ and $\bnabla\times$ stand for the gradient and curl in $\mathbb{R}^{3}$. For a vector field tangent to a plane embedded in $\mathbb{R}^{3}$, this can be simplified to

[TABLE]

where $\bm{N}$ is the normal to the plane and $\Psi$ is a scalar function. Therefore, for a plane embedded in $\mathbb{R}^{3}$, a tangent vector field can be represented in terms of two scalar fields, $\Phi$ and $\Psi$. This property can be generalized to arbitrary surfaces in terms of their intrinsic differential geometry, i.e. not relying on their embedding in $\mathbb{R}^{3}$, as a special case of the Hodge decomposition for $n$-forms (Do Carmo, 1992). A vector field $\bm{v}$ tangent to a surface $\Gamma$ can be decomposed as

[TABLE]

where $\alpha$ and $\beta$ are scalar fields on $\Gamma$ and $\bm{h}$ is a harmonic vector field, satisfying $\bnabla\bcdot\bm{h}=0$ and $\bnabla\times\bm{h}=0$. We note that the curl operator $\bnabla\times$ on a surface, an instance of exterior derivative, is defined differently to its counterpart in Euclidean space. For instance, applied on a scalar function $\bnabla\times\beta$ is a vector with components $(\bnabla\times\beta)^{a}=\epsilon^{ab}\bnabla_{b}\beta$, where $\bm{\epsilon}$ is the antisymmetric tensor

[TABLE]

with $\varepsilon^{ab}$ the Levi-Civita symbols defined by the matrix

[TABLE]

For simply connected surfaces, i.e. closed surfaces with genus equal to 0, there is only a trivial harmonic vector field, $\bm{h}=\bm{0}$, and $\bm{v}$ can be described in terms of the two scalar fields $\alpha$ and $\beta$ (see figure 5 for an example on an ellipsoid).

In the absence of shape changes, from equation (47) it is clear that an inextensible flow satisfies $\bnabla\bcdot\bm{v}=0$. In this case, $\bm{v}$ can be represented in terms of a stream function as $\bm{v}=\bnabla\times\beta$. This approach was introduced by Secomb & Skalak (1982) to describe flows in fluid surfaces with fixed shape and used more recently by various authors (Morris & Turner, 2015; Sigurdsson & Atzberger, 2016; Reuther & Voigt, 2016; Gross & Atzberger, 2018; Mickelin et al., 2018). However, we note that for inextensible surfaces that change shape, both $\alpha$ and $\beta$ need be considered. In this case, using the fact that $\bnabla\bcdot\bm{v}=\bnabla\bcdot\bnabla\alpha+\bnabla\bcdot\bnabla\times\beta=\bnabla\bcdot\bnabla\alpha=\Delta\alpha$, where $\Delta\alpha=\bnabla\bcdot\bnabla\alpha$, it follows from equation (47) that in an inextensible flow $\alpha$ and $v_{n}$ satisfy the constraint

[TABLE]

3 Physical models of fluid surfaces

In this section we examine classical models for fluid deformable surfaces, two used to model lipid bilayers and one applicable to the cell cortex. Thanks to the tools introduced above and Onsager’s formalism, we derive the corresponding governing equations in their full three-dimensional and nonlinear generality.

3.1 Lipid bilayers: An inextensible viscous layer with bending energy

Lipid membranes are interfacial viscous fluids with bending elasticity. The interplay between viscosity and elasticity determines their relaxation dynamics after they are brought out-of-equilibrium by external forces or biological activity. These two essential mechanical features of lipid membranes, their out-of-plane elasticity and interfacial viscosity, have often been examined separately. The mechanical equilibrium of lipid bilayers can be understood to a large extent with the classical bending model of Helfrich (Helfrich, 1973; Lipowsky, 1991; Jülicher & Lipowsky, 1993). For that reason, studies of lipid bilayers at scales beyond tens of nanometers have mainly focused on this model, e.g. in investigations of equilibrium configurations of closed vesicles under geometric constraints, such as fixed surface area or fixed enclosed volume (Steigmann, 1999; Capovilla & Guven, 2002; Tu & Ou-Yang, 2004; Feng & Klug, 2006; Rangarajan & Gao, 2015; Sauer et al., 2017). Beyond the Helfrich model, and subsequent refinements such as the Area Difference Elasticity model (Seifert, 1997), more general models are required to describe the dynamical transformations that bilayers undergo, which should capture the interfacial dissipative mechanisms that dominate at sub-cellular scales. The interfacial hydrodynamics of bilayer membranes was first examined separately from membrane deformation, i.e. assuming fixed membrane shape. These studies focused on the mobility of membrane inclusions, such as proteins, starting with the seminal work of Saffman & Delbrück (1975) on planar lipid bilayers. Subsequent studies have considered the effect of fluid boundaries (Stone & Ajdari, 1998) or the (fixed) shape of the fluid membrane (Levine et al., 2004; Henle & Levine, 2010; Sigurdsson & Atzberger, 2016). Interfacial flows of vesicles induced by shear bulk flows were also considered at fixed vesicle shape (Secomb & Skalak, 1982). Following the seminal works of Scriven (1960) and Aris (1962) on the hydrodynamics of insoluble fluid films, Barthes-Biesel & Sgaier (1985) examined the interfacial flow of vesicles in a shear flow allowing for infinitesimal shape deformations. More recently, a geometrically non-linear model for an inextensible viscous interfacial fluid with bending rigidity was examined, formulated geometrically, and exercised under the assumption of axisymmetry (Arroyo & Desimone, 2009; Arroyo et al., 2010). Along these lines, there is an increasing interest in the community of applied and computational mathematics to develop numerical methods to solve the three-dimensional equations governing inextensible viscous interfaces with curvature elasticity (Nitschke et al., 2012; Rodrigues et al., 2015; Reuther & Voigt, 2016; Barrett et al., 2016b). This model provides a first approximation to the dynamical behaviour of lipid membranes.

Here, we formulate this model based on Onsager’s variational formalism (Arroyo & Desimone, 2009), and derive the Euler-Lagrange equations. We first introduce the bending energy of the bilayer, or Helfrich energy,

[TABLE]

where $\kappa$ and $\bar{\kappa}$ are the bending and Gaussian bending modulus respectively, which we assume to be homogeneous on $\Gamma_{t}$, $H_{0}$ is the spontaneous curvature of the membrane, and $K$ is the Gaussian curvature $K=\det(k\indices{{}^{a}_{b}})$. $\kappa$ is a positive number and $\bar{\kappa}\approx-\kappa$. Thus, this energy penalizes deviations of the mean curvature away from the spontaneous curvature $H_{0}$ and disfavors regions with negative Gaussian curvature (see figure 6). Recalling the Gauss-Bonnet theorem (do Carmo, 2016), according to which $\int_{\Gamma_{t}}KdS$ is a topological invariant, we can ignore the second term in the Helfrich energy for closed surfaces of fixed topology. For simplicity, we restrict our attention to symmetric bilayers, with the same composition in both monolayers, for which $H_{0}=0$. Thus, we can rewrite the Helfrich energy as

[TABLE]

The free energy depends on the set of state variables of the system, usually denoted by $X$. In this case, the material parametrization is the state variable of the system, $X=\{\bm{\phi}\}$ . We note, however, that since the energy only depends on the shape of $\Gamma_{t}$, $\bm{\phi}$ can be replaced by any other ALE parametrization $\bm{\psi}$. In Onsager’s variational principle, dissipative mechanisms are introduced through dissipation potentials. The Newtonian shear rheology of lipid membranes (Dimova et al., 2006) is encoded in the following dissipation potential

[TABLE]

where $\mu$ is the (in-plane) shear viscosity of the monolayer and $|\bm{d}|^{2}=d_{ab}g^{ac}g^{bd}d_{cd}$. Since we assume that the deformation of the membrane is inextensible, only the deviatoric part of $\bm{d}$ matters in the above definition. The dissipation potential depends on the state variable of the system, $\bm{\phi}$, through shape in $dS$ and $\bm{g}$, but primarily on the variable representing the rate of change of the state, $\bm{V}$. The variables that represent the processes that change the state of the system and produce dissipation are called process variables and denoted by $V$; here $V=\{\bm{V}\}$. Here the process variable is simply the time-derivative of the state variable $\bm{V}=\partial_{t}\bm{\phi}$, but this is not necessarily the case. For instance, if we had used $\bm{\psi}$ rather than $\bm{\phi}$ as state variable, then $\bm{W}=\partial_{t}\bm{\psi}$ would not be a meaningful variable to encode dissipation since it does not represent a physical velocity. For thermodynamical consistency, the dissipation potential $\mathcal{D}(X;V)$ must be a convex function of $V$ with minimum at $V=0$ (Arroyo et al., 2018). We further assume that $\mathcal{D}(X;0)=0$, so that $\mathcal{D}(X;V)\geq 0$. If external forces $\bm{F}$ are applied, these introduce a power input

[TABLE]

One can also include the dissipation potential associated to the bulk viscous fluid where the membrane is embedded. Here, we ignore bulk hydrodynamical forces to focus on the fluid membrane, an assumption which is physically justified for phenomena below the Saffman-Delbrück lengthscale $l_{\text{SD}}=\mu/\mu_{b}$, where $\mu_{b}$ is the bulk viscosity (Saffman & Delbrück, 1975; Arroyo & Desimone, 2009). For a lipid membrane, $l_{\text{SD}}\approx$5\text{,}\mathrm{\SIUnitSymbolMicro m}$$.

Onsager’s variational principle establishes a competition between energy release rate, power and dissipation through the Rayleighian, which takes the form

[TABLE]

Here, the rate of change of the energy $D_{t}\mathcal{F}_{H}[\bm{\phi};\bm{V}]$ is

[TABLE]

where we have used that (Capovilla & Guven, 2002; Arroyo & Desimone, 2009)

[TABLE]

Then, Onsager’s principle states that process variables minimize the Rayleighian

[TABLE]

subject to constraints $Q[\bm{\phi};\bm{V}]$. Here, we consider that the surface is inextensible

[TABLE]

Furthermore, due to osmotic effects, it can often be assumed that cells and vesicles maintain their volume constant, and hence

[TABLE]

To enforce these constraints, we can introduce the Lagrangian

[TABLE]

$P$ is the pressure in the vesicle and $\gamma$ is a component of the surface tension. Then, Onsager’s principle subject to constraints can be written as a saddle-point problem

[TABLE]

From the stationarity conditions of equation (67), one finds the weak and the strong form of the governing equations (see appendix D), the latter of which take the form

[TABLE]

Here, $\bm{\Sigma}^{a}$ is the so-called surface stress vector,

[TABLE]

where $\bm{\sigma}$ is the in-plane stress

[TABLE]

and $\bm{\sigma}_{n}$ is a vector of normal stresses

[TABLE]

When multiplied by a unit vector $\bm{l}$ in $T_{\bm{x}}\Gamma_{t}$, $\bm{\Sigma}^{b}l_{b}$ is the three-dimensional force per unit length across a curve passing through $\bm{x}$ and perpendicular to $\bm{l}$. Note that $\sigma_{n}^{a}$ represents bending moments caused by curvature imbalances. Finally, $\bm{B}$ is the field of (external) body forces

[TABLE]

Eqs. (68)-(70) express balance of linear momentum and conservation of mass (inextensibility) and enclosed volume in a fully nonlinear regime. Alternatively, one could have derived equation (68) from local force balance on the membrane as in (Salbreux & Jülicher, 2017), and postulated the constitutive laws Eqs. (72) and (73). Thus, by starting from different ingredients (a Rayleighian expressing energy release-rate, dissipation and power input) and invoking a variational principle subject to constraints, Onsager’s formalism recovers these equations in a systematic and transparent way. As a direct corollary of Onsager’s principle, it is easy to see that, in the absence of external power inputs, the free energy is a Lyapunov functional of the dynamics, i.e. $\mathcal{F}$ is a decreasing functional (Arroyo et al., 2018)

[TABLE]

which provides a nonlinear notion of stability for the dynamics. From a computational point of view, only the weak form of the stationarity conditions issuing from Onsager’s formalism is required for a space discretization based on finite elements. Onsager’s formalism also provides a framework to formulate nonlinearly stable variational time-integrators, as described in section 4.

We finally note that the choice of state and process variables is not unique. For instance, as mentioned earlier, we can choose an ALE parametrization $\bm{\psi}$ instead of $\bm{\phi}$ as state variable, since the free energy depends only on shape. The velocity $\bm{V}$, our process variable, then needs to be split into a normal and a tangential components $\bm{V}=\bm{v}+v_{n}\bm{N}$, where $v_{n}=\bm{W}\bcdot\bm{N}$. More specifically and for the ALE parametrization in equation (48), we have $v_{n}=\partial_{t}h\bm{M}\bcdot\bm{N}$. We can then rewrite the Lagrangian in terms of $\partial_{t}h$ and $\bm{v}$. We can further decompose $\bm{v}=\bnabla\alpha+\bnabla\times\beta$. The governing equations issued from any of these choices look very different but describe the same dynamics. For instance, the dynamics obtained from

[TABLE]

are equivalent to those resulting from equation (67). While the choice of variables $X=\{\bm{\phi}\}$ and $V=\{\bm{V}\}$ is natural from a modelling viewpoint, the choice $X=\{h\}$ and $V=\{v_{n},\alpha,\beta\}$ is better suited from a computational viewpoint, as will become clear in section 4.

3.2 Lipid bilayers: The Seifert-Langer model

The previous model provides a first approach to the mechanics of lipid bilayers. It is often overlooked, however, that by ignoring the bilayer architecture it fails to capture many important phenomena. Seifert and Langer developed a continuum model explicitly accounting for the bilayer architecture and capturing the major energetic driving forces and dissipative drag forces involved in the dynamics of lipid membranes (Seifert & Langer, 1993). The elastic forces in this theory appear in response to bending of the membrane, as in the previous model, but also to monolayer stretching (see figure 7). As viscous effects, the in-plane Newtonian rheology of the lipid bilayer (Dimova et al., 2006) is included through shear and dilatation dissipations, and the frictional coupling between the two monolayers opposing inter-monolayer slippage is also included. This model provided predictions about the relaxation dynamics of membrane fluctuations. Importantly, its material parameters can be experimentally measured (Dimova et al., 2006). The work of Seifert & Langer (1993), along with (Evans & Yeung, 1994), highlighted the role of inter-monolayer friction as a “hidden” but significant dissipative effect. This physical model was originally introduced and has been predominantly exercised under the restricted assumptions of linearized disturbances around a planar state (Seifert & Langer, 1993; Fournier, 2015) or of simplified and fixed membrane shape (Evans & Yeung, 1994). These approximations, however, hide much of the interaction between shape dynamics and interfacial hydrodynamics, which is mediated by membrane curvature. This was demonstrated by the linearization of the theory about spherical or cylindrical configurations (Rahimi, 2013) and by simulations based on a fully non-linear version of this theory, albeit axisymmetric (Rahimi & Arroyo, 2012), which further demonstrated the geometry-dependent subtle interplay between all the ingredients in figure 7 at multiple scales. Seifert and Langer’s (SL) model is conceptually simple, captures sufficient physics to describe a plethora of dynamical phenomena, and can be the basis for more sophisticated dynamical models including for instance lipid tilt near molecular inclusions (Hamm & Kozlov, 2000, 1998) or the physicochemical interaction of lipids with scaffolding or integral proteins (Brochard-Wyart & de Gennes, 2002; Arroyo et al., 2018). Here we formulate and develop numerical calculations with this model in a three-dimensional and fully non-linear setup which, to our best knowledge, has not been examined before.

In this model, $\Gamma_{t}$ characterizes the bilayer mid-surface (see figure 8). In addition to the Helfrich energy of the form of equation (57), the Seiftert-Langer model accounts for the stretching elasticity of each of the monolayers through the functional

[TABLE]

where the fields $\breve{\rho}^{\pm}$ represent the lipid density at the neutral surface of the upper (+) and lower (-) monolayers, measured in units of the equilibrium density, which differ from the lipid density at the bilayer mid-surface according to

[TABLE]

with $d\approx$1\text{,}\mathrm{nm}$$ the distance between the neutral surfaces and the bilayer mid-surface. Unless otherwise noted, a functional containing $\pm$ implies a summation on the $+$ and $-$ monolayers. For convenience, in this section we use an Eulerian parametrization $\bm{\chi}$ to derive the governing equations. The free energy depends on $\bm{\chi}$ and the two density fields $\rho^{+}$ and $\rho^{-}$, representing the density of lipids on the mid-surface, and thus $X=\{\bm{\chi},\rho^{+},\rho^{-}\}$. We take into account three main dissipation mechanisms in the bilayer. First we consider the dissipation due to in-plane shear in each monolayer, which takes the form

[TABLE]

where $\mu$ is the shear viscosity and

[TABLE]

is the rate-of-deformation tensor for each monolayer (see equation (30)). We consider three process variables, $v_{n}=\partial_{t}\bm{\chi}\bcdot\bm{N}$, which determines shape changes, and $\bm{v}^{+}$ and $\bm{v}^{-}$, which determine the tangential flow of lipids in each monolayer. Thus, $V=\{v_{n},\bm{v}^{+},\bm{v}^{-}\}$. Additionally, we consider a dilatational dissipation

[TABLE]

where $\lambda$ is the dilatational viscosity. Finally, we consider the inter-monolayer friction caused by the relative slippage of one monolayer with respect to the other

[TABLE]

where $b_{I}$ is the inter-monolayer friction coefficient. Thus, the total dissipation is

[TABLE]

The rate of change of the free energy is

[TABLE]

Note carefully that $D_{t}\mathcal{F}_{S}\left[\bm{\chi};\rho^{\pm};v_{n},\partial_{t}\rho^{\pm}\right]$ depends on $\partial_{t}\rho^{\pm}$ rather than on the process variables $\bm{v}^{\pm}$. We invoke the equations encoding conservation of mass

[TABLE]

where $\Pi$ is referred to as a process operator, to express $D_{t}\mathcal{F}_{S}$ in terms of the process variables, in equal footing with $\mathcal{D}$, towards applying Onsager’s formalism (Rahimi & Arroyo, 2012). Process operators, usually linear operators, relate the rate of change of state variables, in this case ${\partial}_{t}\rho^{\pm}$, with process variables, $v_{n}$ and $\bm{v}^{\pm}$. In general, we write

[TABLE]

In the previous model, the process operator was trivial $\dot{X}=V$ ($\partial_{t}\bm{\phi}=\bm{V}$). Using equation (85), the rate of change of the energy can be written as

[TABLE]

and the Lagrangian

[TABLE]

where here

[TABLE]

Then, Onsager’s variational principle states that

[TABLE]

The stationarity conditions issued from Onsager’s principle provide equations for $P$ and the fields $v_{n}$ and $\bm{v}^{\pm}$. To find the time-evolution of the density fields $\rho^{\pm}$, the process operator (equation (85)) needs to be integrated in time. We stress that Onsager’s variational principle provides directly the weak form of the problem, which can be directly discretized with finite elements. For completeness, we derive using Onsager’s formalism the stress tensor and strong form of the governing equations of SL model, which to the best of our knowledge have not been presented before in the fully nonlinear case. The tangential and normal components of the stress of each monolayer can be identified as (see appendix E)

[TABLE]

and

[TABLE]

Note that, aside from the terms involving $k_{S}$, the expressions are similar to those of previous model. Density imbalances generate a source of in-plane stress, but also lead to bending moments. Funthermore, the bending rigidity of the bilayer is $\kappa+2k_{S}(\rho^{\pm}d)^{2}$, which includes the effect of Helfrich and stretching energies. Balance of linear momentum tangent to the surface on the upper and lower monolayers reads

[TABLE]

where $b_{I}\left(\bm{v}^{+}-\bm{v}^{-}\right)$ identifies the force exerted by the lower monolayer on the upper monolayer due to intermonolayer friction. Finally, balance of linear momentum perpendicular to the bilayer leads to

[TABLE]

Seifert and Langer first introduced a linearized version of these equations around a planar state (Seifert & Langer, 1993), which has been recently reviewed in the context of Onsager’s principle (Fournier, 2015). The stress tensors in equation (91) and equation (92) are similar to those found in (Rahimi & Arroyo, 2012), using the Doyle-Ericksen formula of continuum mechanics. Our general and systematic derivation shows the ability of Onsager’s formalism to derive complex models mixing different physics in a fully non-linear setting, which would otherwise be difficult to rationalize. For instance, although not unconceivable, it is difficult to postulate the constitutive relation for the in-plane stress in equation (91).

3.3 The cell cortex: A viscous layer driven by active tension

The cell cortex is a layer of cross-linked actin filaments lying just beneath the plasma membrane of animal cells (Bray & White, 1988). The thickness of this layer is of hundreds of nanometers, while the typical size of an animal cell is of tens of microns. Thus, this layer can be considered as a quasi two-dimensional material. In addition to actin, this network is crowded with polymerization regulators, cross-linkers, or myosin motors, which bind to actin filaments. By consuming ATP, these molecular motors pull on actin filaments and generate active tension. In turn, this active tension, if non-uniform, generates actin flows and drives shape changes (Salbreux et al., 2012). Another important property of this actin network is that undergoes dynamic remodelling, with a continuous turnover by polymerization and depolymerization of actin filaments and binding and unbinding of cross-linking proteins (Howard, 2001). This process is characterized by a time-scale in the order of a few tens of seconds. At time-scales shorter than the turnover time, the cortex behaves as an elastic network. At longer time-scales, the dynamic remodelling of the cortex leads to a fluid-like viscous behaviour with active tension.

Following previous works (Turlier et al., 2014; Prost et al., 2015), we consider an active gel model of the cortex as an isotropic viscous material with active tension confined to a surface and undergoing turnover, in which viscosity, active tension, and depolymerization depend on the thickness of the cortex. This model can describe phenomena at time-scales of a minute and longer, where elastic energy storage in the network becomes negligible, and does not account for the architecture of the network, e.g. the orientation of the actin filaments, which may not be appropriate in some important examples such as during cytokinesis (Reymann et al., 2016). Furthermore, we assume that the viscous forces exerted by the cytosol and the external fluid medium are negligible. Using common estimates of cortex 2D viscosity ($\mu=27\times 10^{-4}$\text{,}\mathrm{Pa}\text{,}\mathrm{s}\text{,}\mathrm{m} (Bergert *et al.*, [2015](#bib.bib13))), an estimate for the Saffman-Delbrück length scale is $l_{\text{DS}}\approx$3\text{\,}\mathrm{m}. Thus, and given the size of cells, neglecting bulk viscosity is well-justified. Turlier et al. (2014) and previous works (Bergert et al., 2015; Saha et al., 2016) were restricted to axisymmetric or to two-dimensional configurations, and derived the active gel equations from the stress tensor and force balance. Here, we develop a fully three-dimensional and geometrically non-linear version of this active gel model, and derive the governing equations using Onsager’s formalism.

Mathematically, we characterize the cortex as a fluid surface $\Gamma_{t}$, described here for the purpose of deriving the governing equations with a Lagrangian parametrization $\bm{\phi}$, with a space-varying thickness $\rho$, see figure 9. $\bm{\phi}$ and $\rho$ are our state variables. The process variable in this problem is the velocity field $\bm{V}$ of actin, with a tangential component $\bm{v}$, characterizing the flow of actin on $\Gamma_{t}$, and a normal component $v_{n}\bm{N}$ describing the change of shape of the actin cortex. The viscous rheology of the cortex is characterized by a dissipation potential, similar to that of lipid bilayers

[TABLE]

where here $\mu$ is the bulk shear viscosity of the cortex. This dissipation potential can be obtained by integrating over the thickness a three-dimensional shear dissipation potential $\int\mu|\bm{D}|^{2}dV$, with $\bm{D}$ the three-dimensional rate-of-deformation tensor, for an incompressible slab of gel with a plain stress assumption, i.e. assuming that $\bm{D}=d_{ab}\bm{e}^{a}\otimes\bm{e}^{b}+D_{nn}\bm{n}\otimes\bm{n}$ and $D_{nn}=-\text{tr}\bm{d}$ (Salbreux et al., 2009; Turlier et al., 2014). To introduce the active tension generated by the activity of myosin motors, we consider a power input of the form

[TABLE]

where $\xi$ is a measure of myosin activity, which may depend on cortical density $\rho$, see discussion in section 5.3. This leads to an active surface tension $\gamma=\xi(\rho)\rho$. Since $\text{tr}\bm{d}$ measures the rate at which local area expands (positive $\text{tr}\bm{d}$) or contracts (negative $\text{tr}\bm{d}$), for a positive $\gamma$ the power input functional will drive the contraction of cortex area. As we neglect the elastic behaviour of the cortex, there is no free energy associated to the problem. Introducing a cell volume constraint, we obtain the Lagrangian

[TABLE]

and the dynamics follows from

[TABLE]

From the Euler-Lagrange equations we identify the constitutive law

[TABLE]

and the statement of balance of linear momentum, this time in the absence of bending moments,

[TABLE]

with the last equation generalizing Laplace’s law. To relate the rate of change of $\rho$ and $\bm{V}$, we consider balance of cortex material

[TABLE]

where the first term in the right hand side stands for actin polymerization, which, since polymerization nucleators are located at the plasma membrane, is assumed to occur at a constant rate $k_{p}$ independent of the thickness, and the second term stands for actin depolymerization, which is proportional to the local thickness, $k_{d}\rho$. The ratio $\rho_{0}=k_{p}/k_{d}$ determines the thickness at steady-state. By defining the characteristic turnover time as $\tau=1/k_{d}$, we can rewrite the previous equation as

[TABLE]

4 Discretization of the mechanics of fluid surfaces

In this section, we introduce a general discretization framework for the simulation of fluid surfaces. First, we introduce a variational time-integrator based on Onsager’s principle, which is unconditionally stable by construction. Then, we introduce the spatial discretization of the different fields on the time-evolving surface. We end with the derivation of the discrete equations for an inextensible fluid surface with bending elasticity as a reference example.

4.1 Time discretization: Variational time-integrator based on Onsager’s principle

To integrate in time the dynamics of continuum mechanical systems, a common approach is to first discretize in space, obtain a system of ordinary differential equations, which is then integrated in time with specialized algorithms. The fact that the dynamics in the models examined here emerge from a variational principle provides an alternative approach: to discretize in time the variational principle itself. Time-integrators based on the discretization of a variational principle are usually referred to as variational time-integrators, and have been widely employed, for instance, for the discretization of Hamilton’s principle in conservative systems including molecular dynamics (Frenkel & Smit, 2001) and elastodynamics (Lew et al., 2004), and in the context of dissipative systems (Ortiz & Stainier, 1999; Peco et al., 2013). Variational time-integrators inherit qualitative properties of the associated time-continuous problem. For instance, in the case of time-integrators based on Hamilton’s principle, Noether’s theorem ensures that symmetries in the discrete action result in conserved currents as in the original continuous theory. Here, we propose a first order variational time-integrator for Onsager’s principle that inherits that $\mathcal{F}$ is a Lyapunov functional of the dynamics, see equation (75). This feature provides nonlinear stability to the resulting discrete dynamics by construction.

We consider here a general statement of Onsager’s variational principle, with a set of state variables $X$, a set of process variables $V$, obeying Onsager’s principle in equation (63) and a process operator as in equation (86). For simplicity, we neglect constraints in our discussion but they can be added by substituting the Rayleighian by the corresponding Lagrangian without changing the essence of the proposed variational integrator. Let us consider a time discretization $\{t^{1},\dots,t^{N}\}$ and let us start with a trivial process operator $\partial_{t}{X}=V$. We will consider here the simplest low order version of implicit variational time-integrator based on Onsager’s principle, and leave the investigation of higher-order schemes to future work. We approximate $V^{n}=V(t^{n})$ with a simple backward difference

[TABLE]

where $X^{n}=X(t^{n})$ and $\Delta t^{n}=t^{n+1}-t^{n}$. The dissipation potential and the power can now be approximated as

[TABLE]

To discretize the Rayleighian, we also need to discretize the rate of change of the free energy. Rather than resorting to an expression like $\dot{\mathcal{F}}\approx D\mathcal{F}\bcdot\left(X^{n+1}-X^{n}\right)/\Delta t^{n}$, we consider

[TABLE]

or a similar higher-order finite difference. This approach ensures that $\mathcal{F}$ is a Lyapunov functional of the dynamics, as we prove below, and retains the full non-linearity of $\mathcal{F}$ in the formulation. Using the previous expressions we define the discrete Rayleighian as

[TABLE]

where we have ignored the constant term $\mathcal{F}(X^{n})/\Delta t^{n}$. Then, the incremental Onsager’s principle is given by

[TABLE]

Thus, the dynamical problem arising from our variational time-discretization can be interpreted as an energy minimization problem for $\mathcal{F}$, which is usually a non-linear function of $X^{n+1}$, with the addition of a convex (and often quadratic) function of $X^{n+1}$, $\mathcal{D}$, subject to the external forces represented in $\mathcal{P}$. The weight of $\mathcal{F}$ relative to $\mathcal{D}$ is controlled by $\Delta t^{n}$, which can be decreased to ease the solvability of the problem by increasing the influence of the convex functional $\mathcal{D}$, or increased to allow the system to reach equilibrium faster. Let us now prove that, for a homogeneous problem ($\mathcal{P}(X;V)=0$), the free energy is a Lyapunov functional of the dynamics. We evaluate the Rayleighians

[TABLE]

where we have used that $\mathcal{D}\left(X^{n};0\right)=0$, as discussed in section 3.1. Since $X^{n+1}$ minimizes $\mathcal{R}^{n}$, it is clear that $\mathcal{R}^{n}\left(X^{n};X^{n+1}\right)-\mathcal{R}^{n}\left(X^{n};X^{n}\right)\leq 0$. Then,

[TABLE]

where we have used that $\mathcal{D}\left(X^{n};\frac{X^{n+1}-X^{n}}{\Delta t^{n}}\right)$ is positive. Therefore, we obtain

[TABLE]

which shows that $\mathcal{F}$ is a Lyapunov functional of the discrete dynamics. Thus, the time-step is not limited by stability, but rather by accuracy and solvability of the non-linear optimization problem in equation (107), which becomes “easier” or “more convex” for small $\Delta t^{n}$. The ability to take stably large time-steps is particularly useful in stiff problems, such as those involving the Helfrich curvature energy.

When the process operator is not trivial, i.e. $\partial_{t}{X}\neq V$, the approach above needs to be modified. For those cases, we can keep $V^{n+1}$ as the variable of the discrete Onsager’s principle and discretize the process operator in different ways. As a first approach, we can consider a simple forward Euler approximation for the process operator

[TABLE]

We can then rewrite equation (105) as

[TABLE]

This approximation still retains the non-linearity of $\mathcal{F}$ and is thus implicit in this sense. We can now define the Rayleghian as

[TABLE]

and solve

[TABLE]

Finally, we can recover $X^{n+1}$ from equation (111). With this simple forward approximation for the process operator, however, the accuracy and stability of the integration can be very limited. As a better alternative, we consider a backward Euler approximation of the process operator, which involves solving

[TABLE]

together with the minimization of the Rayleighian

[TABLE]

That is, one needs to solve the system

[TABLE]

for $X^{n+1}$ and $V^{n+1}$ simultaneously. It is easily shown that with any of these discretizations, $\mathcal{F}$ is also a Lyapunov function of the dynamics in the absence of power input, thus retaining the nonlinear stability of the time-discretization scheme.

4.2 Spatial discretization

In this section, we examine the spatial discretization of $\Gamma_{t}$ and the different fields defined on it. For simplicity, let us start by examining the numerical parametrization of a generic surface $\Gamma$. We first note that, since models for fluid surfaces usually involve the shape operator $\bm{k}$, this tensor needs to be square-integrable on $\Gamma$. For that reason, the parametrization of $\Gamma$ must be a square-integrable function with square-integrable first- and second-order derivatives; we call such a surface a $H^{2}$ surface. The problem of discretizing a $H^{2}$ surface may be addressed resorting to different numerical frameworks, such as higher-order B-splines as in isogeometric methods (Piegl & Tiller, 2012; Sauer et al., 2017) or max-ent approximants (Millán et al., 2011). Another versatile technique to discretize smooth surfaces based on meshes with arbitrary connectivity is subdivision surfaces. Here we focus on Loop subdivision surfaces based on triangular meshes (Loop, 1987; Stam, 1999; Biermann et al., 2000; Cirak et al., 2000; Cirak & Ortiz, 2001; Cirak & Long, 2011; Torres-Sánchez, 2017). To define the discretization of $\Gamma$ with subdivision surfaces, we consider a control mesh made of triangles $E=1,\dots,N_{e}$ whose edges join the set of control points with positions $\{\bm{x}_{I}\}_{I=1}^{N_{n}}$. For each triangle in the mesh, we define the parametrization $\bm{\psi}^{E}(\bm{\xi}):\tilde{\Gamma}\rightarrow\mathbb{R}^{3}$, with $\tilde{\Gamma}$ the reference triangle (see figure 10), by

[TABLE]

where $B_{I}^{E}$ represents the subdivision basis function associated to node $I$ at element $E$ and $\langle E\rangle^{1}$ identifies the first ring of nodes surrounding the element, including the nodes forming the element and all first neighbours to them. We denote by $\Gamma^{E}=\bm{\psi}^{E}\left(\tilde{\Gamma}\right)$ the curved triangle obtained by the local parametrization in equation (118). It can be shown that these curved triangles are disjoint (except at the edges) and that their union $\Gamma=\cup_{E=1}^{N_{e}}\Gamma^{E}$ defines a $C^{2}$-continuous surface almost everywhere, except at a finite number of points where it is $C^{1}$. These points coincide with the image of irregular nodes in the control mesh, which are those with a connectivity different from 6. There, the surface is continuous with continuous derivative but presents a discontinuity in the second derivative. Thus, from the perspective of differential geometry, the set $\cup_{E=1}^{N_{e}}\left(E,\tilde{\Gamma},\psi^{E}\right)$ defines an atlas of charts that parametrize the $H^{2}$ surface $\Gamma$.

Now we consider ALE parametrizations of the form derived in section 2.5.1. For the parametrization of the surface $\Gamma_{t_{0}}$, we write

[TABLE]

Thus, the control mesh in our scheme is given by the position of the control points $\left\{\bm{x}_{0I}\right\}_{I=1}^{N_{n}}$. We also define the fields

[TABLE]

The parametrization of the deformed surface $\Gamma_{t}$ then reads

[TABLE]

We note that, given that $\bm{\psi}_{0}$, $h$ and $\bm{M}$ are in $H^{2}$, $\bm{\psi}$ is also in $H^{2}$. We also note that, if we had used the normal to the reference surface $\bm{N}_{0}$ instead of $\bm{M}$, because the calculation of $\bm{N}_{0}$ already involves first order derivatives of $\bm{\psi}_{0}$, we would need $\Gamma_{t_{0}}$ be $C^{2}$ everywhere, which cannot be achieved with subdivision surfaces. This is the reason why we choose the field of directors as in equation (121), where $\bm{M}_{I}$ can be chosen to approximate the true field of normals, for instance in a least-squares sense.

We can consider other kinds of basis functions. In particular, we consider the set of linear basis functions $N_{I}(\bm{\xi})$ with $I\in\langle E\rangle^{0}$, the zeroth-ring of nodes of the element, defined by

[TABLE]

where $E_{1},E_{2}$ and $E_{3}$ denote the labels of the three nodes forming the element $E$. We can then discretize fields on $\Gamma_{t}$ with $N_{I}$ if they only need to be in $H^{1}$ (that is, square-integrable functions with a square-integrable derivative). For instance, a density field, which only appears in the free energy and dissipation potentials through its value and first order derivatives, can be discretized as

[TABLE]

For notational simplicity in this and following equations, we write $N_{I}^{E}$ where we should write $N_{I}^{E}\circ\left(\bm{\psi}^{E}\right)^{-1}$ since $N_{I}^{E}$ take values in the parametric domain $\tilde{\Gamma}$ whereas $\rho$ is a field on $\Gamma_{t}$. Following the same arguments, we could be tempted to discretize the components of $\bm{v}$ as

[TABLE]

since the Rayleighian also depends on $\bm{v}$ through its value and its first-order derivatives only. This, however, requires that a basis $\bm{e}_{a}$, continuous and with first square-integrable derivatives, is defined everywhere on the surface so that $\bm{v}=v^{a}\bm{e}_{a}$ is continuous and with square-integrable derivatives. However, one cannot define such a basis for a closed surface as a consequence of the hairy ball theorem (for instance, polar coordinates in the sphere present singularities at the poles). Using the canonical basis of the parametrization, we could try to discretize $\bm{v}=v^{a}\partial_{a}\bm{\psi}^{E},$ but $\partial_{a}\bm{\psi}^{E}$ is discontinuous across elements due to the jump in the definition of local coordinates. A possible solution to this problem is to increase the number of degrees of freedom used to describe $\bm{v}$ and discretize the three components of $\bm{v}$ in the global basis of Euclidean space

[TABLE]

Being the basis vectors $\bm{i}_{\alpha}$ constant, we could discretize

[TABLE]

However, $\bm{v}$ being tangent to $\Gamma_{t}$, its three components $v^{1},v^{2}$ and $v^{3}$ are not independent and one would need to introduce the additional constraint $\bm{v}^{E}\bcdot\bm{N}=0$ such as in (Fries, 2018; Reuther & Voigt, 2018). A more convenient option is to recall the Hodge decomposition of $\bm{v}$ in equation (52) and discretize the scalar fields $\alpha$ and $\beta$. We note that $\alpha$ and $\beta$ need to be in $H^{2}$ for $\bm{d}$ to be well-defined, and for this reason we use subdivision basis functions to discretize them

[TABLE]

Apart from $h$ and the vector potentials $\alpha$ and $\beta$, in some models we need to discretize Lagrange multiplier fields, such as the surface tension $\gamma$. Since $\gamma$ acts as a Lagrange multiplier, the space of basis functions for $\gamma$ needs to be chosen with care to ensure that the discretization satisfies the discrete inf-sup condition (Brezzi & Fortin, 2012).

Similarly to previous works in isogeometric analysis (Dortdivanlioglu et al., 2018), we consider a macro-element approach where Lagrange multipliers are approximated using a coarsened mesh. We consider a coarse mesh of triangles $E^{c}\in 1,\dots,N_{e}^{c}$, and subdivide each triangle $E^{c}$ into four triangles $(4E^{c},4E^{c}+1,4E^{c}+2,4E^{c}+3)$; the nodes of these meshes are $N_{n}^{c}$ and $N_{n}^{f}$ respectively. This subdivision can be regarded as a map between the parametric domains of two atlases $\cup_{E=1}^{N_{e}^{c}}\left(E,\tilde{\Gamma}\right)$ and $\cup_{E=1}^{N_{e}^{f}}\left(E,\tilde{\Gamma}\right)$, with $N_{e}^{f}=4N_{e}^{c}$, defined by

[TABLE]

see figure 11. We note that the function $O$ does not depend on the positions of the nodes of the mesh. The finer atlas $\cup_{E^{f}=1}^{N_{e}^{f}}\left(E^{f},\tilde{\Gamma}\right)$ is then used to discretize the geometry of the surface and the vector potentials, following Eqs. (122), (128) and (129). On the other hand, we discretize Lagrange multiplier fields with linear elements defined by equation (123) in the coarser atlas

[TABLE]

4.3 Finite element formulation of Onsager’s principle

Here, we show the application of our methodology, based on the variational time-integrator described in section 4.1 and on the space discretization described in section 4.2, to the model of an inextensible viscous fluid surface with bending elasticity (section 3.1). We define the following vectors of nodal coefficients

[TABLE]

containing the degrees of freedom describing the offset, the irrotational and solenoidal vector potentials, and the surface tension. The discrete Lagrangian, now a function, can then be written as

[TABLE]

where the explicit form for the different terms can be found in appendix F. The discrete version of Onsager’s variational principle then leads to the saddle-point problem

[TABLE]

the stationarity conditions of which form a non-linear algebraic system of equations, which we solve using Newton’s method.

4.4 Restricting rigid body motion in simulations

The simulation of fluid surfaces lacking of interaction with the surrounding viscous fluid requires of restricting rigid body motions of the interface since these do not dissipate energy or affect the free energy of the system. To restrict these motions, we impose three translational constraints

[TABLE]

and three rotational constraints

[TABLE]

using six additional Lagrange multipliers.

4.5 Mass conservation: Stabilized finite element formulation

We now address the discretization of mass conservation, which is required in the Seifert-Langer model of lipid bilayers as well as for the simulation of the cell cortex. Since we consider an ALE parametrization, we need to discretize the ALE version of equation (46). We consider an implicit backward Euler scheme in time for this advection-reaction equation. For its space discretization, we consider a stream-upwind Petrov Garlerkin (SUPG) method (Donea & Huerta, 2003), which treats the convective term by adding controlled numerical diffusion in a consistent manner. The equations of conservation of mass and balance of linear momentum are solved monolithically using Newton’s method. See appendix G for details.

5 Representative simulations of fluid surfaces

In this section, we revisit the models developed in section 3. We use the numerical framework described in the previous section to simulate these models under different conditions, which exemplify the mechanical behaviour of lipid bilayers and the cell cortex.

5.1 Lipid bilayers: An inextensible viscous layer with bending energy

Example 1: Relaxation dynamics from a non-equilibrium non-axisymmetric shape. We first simulate the behaviour of an inextensible viscous layer with curvature elasticity as a first approach to model the elasto-hydrodynamics of lipid bilayers. To test the performance of the numerical methods described in the previous section, we first examine the relaxation of an out-of-equilibrium (and non-axisymmetric) shape given by

[TABLE]

with $R=$100\text{,}\mathrm{\SIUnitSymbolMicro m}, $\lambda_{1}=0.7$ and $\lambda_{2}=0.3$ (see figure [12](#S5.F12)). Using common estimates for the model parameters (Dimova *et al.*, [2006](#bib.bib30); Rahimi & Arroyo, [2012](#bib.bib90)), we choose $\kappa=10^{-19}$\text{\,}\mathrm{J}, $\mu=10^{-9}$\text{,}\mathrm{J}\text{,}\mathrm{s}\text{,}\mathrm{m}$^{-2}$. As expected for a dissipative system in the absence of external inputs, the free energy $\mathcal{F}$ decreases monotonically with time (figure 12A) by dissipating energy (figure 12B). Note that, because of the semi-logarithmic scale, it is difficult to appreciate in Figs. 12A and B that the negative of the rate of change of free energy is equal to the rate of dissipation. The initial shape given in figure 12I relaxes through different non-equilibrium states, figure 12II and figure 12III, until reaching the final equilibrium shape, figure 12IV. In the left panel of Figs. 12I-IV we show the velocity field, which has been split for visualization purposes into its normal (colormap) and tangential (arrows) components. In the right panel of these figures, we show the Lagrange multiplier field $\gamma$ representing the contribution to surface tension of the inextensibility constraint, which shows a smooth behaviour, suggesting that the macro-element approach described in section 4.2 satisfies the discrete inf-sup condition; a more detailed study of this specific will be presented elsewhere. In figure 13A, we compute the $L_{2}$ norm of $\text{tr}\bm{d}$ for three different levels of refinement, marked in blue (with average triangle side $h/R\approx 1/17$, $N_{e}=8128$ and $N_{n}=4066$), green ($h/R\approx 1/34$, $N_{e}=32512$ and $N_{n}=16258$) and red ($h/R\approx 1/68$, $N_{e}=130048$ and $N_{n}=65026$) as a function of time, which measures the error committed in the enforcement of inextensibility. Initially, we observe that the error converges linearly in a log-log scale as the mesh is refined. Even though we use an ALE method to reduce mesh distortion, the dramatic shape changes during the relaxation dynamics require four full remeshing operations, which are marked with yellow circles in Fig. 13A. The resulting meshes are shown in panel of figure 13I-IV and in Movie 1 for the coarser refinement level. To remesh, we follow a three-step procedure. First, we update the reference surface $\Gamma_{t_{0}}$ to $\Gamma_{t}$ using least-squares. This reference surface serves then as a seed for the remeshing algorithm implemented in the VMTK library (Antiga et al., 2008), which assigns an element area following

[TABLE]

where $A_{0}$ is a reference element area for a planar patch, and $c_{s}$ specifies the sensitivity to curvature (in these simulations, $c_{s}=0.1$). Finally, another least squares fit is performed to parametrize the new surface to the initial one $\Gamma_{t}$, which finally sets the new $\Gamma_{t_{0}}$. We note that the first least-squares fit is only performed to give a seed to the meshing algorithm; the essential least-squares fit, using the initial shape as a seed to fit the geometry of a parametrization based on the new mesh, is performed after remeshing. We observe that remeshing increases the error associated to local inextensibility noticeably, but this error remains small. Thus example illustrates the benefit of the ALE method to reduce the frequency of remeshing events. We finally note that the relative error in total area and volume conservation is smaller than 0.1% over the whole dynamics, see figure 13B. The error in volume conservation is very small ($<10^{-11}\%$) until the first remeshing step, where the error presents a jump. This illustrates the success of our non-linear method to impose volume conservation, see section 4.3. On the other hand, it shows the lack of explicit control on volume (and area) conservation during remeshing, which could be incorporated into the least-squares procedure underlying remeshing. Errors in area conservation are smoother in time and larger in magnitude, since it is imposed weakly in terms of local area conservation based on the discretization of $\text{tr}\bm{d}$ and the Lagrange multiplier $\gamma$.

Example 2: Dynamics following hyper-osmotic shocks. As a second example, we examine the effect of osmotic shocks in vesicles. Cells and vesicles are often exposed to changes in the inner and outer chemical composition, which create flows of water through the semipermeable lipid membrane, increasing or decreasing their enclosed volume, and generating shape changes (Staykova et al., 2013; Kosmalska et al., 2015). Here, we simulate the effect of a hyper-osmotic shock by decreasing the enclosed volume at different deflation rates. We start with the equilibrium shape of the previous example using the finest mesh (figure 14-0), and apply a deflation rate of 10\text{\,}\mathrm{nm}$^{3}$\text{\,}\mathrm{ns}$^{-1}$. In a plot comparing the elastic energy stored during deflation and the total volume decrease (blue curve in figure [14](#S5.F14)), we observe a linear dependence. In fact, at this rather small deflation rate, we observe that the shape of the vesicle (figure [14](#S5.F14)A1-A2) follows a sequence of prolate shapes for the given area and volume that are equivalent to those found at equilibrium (Feng & Klug, [2006](#bib.bib40)). We observe, however, a small fluctuation of normal and tangential velocities in the equator of the vesicle, which are a signature of a non-equilibrium symmetry-breaking process. These deviations from axisymmetry become more noticeable at higher deflation rates. For instance, for a deflation rate of 100\text{,}\mathrm{nm}$^{3}$\text{,}\mathrm{ns}$^{-1}$, we observe that the shape starts to deviate from quasi-equilibrium path and velocity variations disturbing axisymmetry are very pronounced (see figure 14-B1), leading to a very different shape as compared to the equilibrium one for the same volume decrease (see figure 14-B2). In agreement with this, we observe that the energy stored during this faster deflation is now higher (green curve). The viscous dissipation of the lipid membrane becomes increasingly dominant as deflation rate increases (see figure 14-C and D respectively). Similarly, we observe that the final shape gets further away from the equilibrium shape, by storing much more elastic energy for a given amount of volume decrease (red and magenta curves). Mechanically, the viscous dissipative forces can be interpreted as a dynamical confinement for the elastic membrane, causing it to transiently buckle and break symmetry.

5.2 Lipid bilayers: Seifter-Langer model

In this section we examine the response of the Seifert-Langer model to monolayer density imbalances, which may arise from chemical perturbations. Membranes in cells and organelles are often exposed to changes in their local lipid density. For instance, proteins and other membrane inclusions, such as polymers, insert in the membrane and locally change the lipid packing (Shibata et al., 2009; Tsafrir et al., 2003). Chemical signals, such as pH disturbances (Khalifat et al., 2008; Fournier et al., 2009), can also alter lipid packing. Furthermore, changes in the local density can occur asymmetrically, affecting only one of the two monolayers. Local density perturbations lead to transient dynamics, where lipid flows and shape changes are tightly coupled and dictated by the interplay between stretching, bending, shear and intermonolayer friction. Thus, these processes constitute an excellent example of application of our theoretical and computational framework. Furthermore, these processes have been previously examined under the assumption of axisymmetry (Rahimi & Arroyo, 2012), which can be used as a reference to verify our numerical procedure.

Following Rahimi & Arroyo (2012), we examine deflated spheroidal prolate vesicles, initially at equilibrium, to which we apply a density disturbance. To prepare the initial state, we start with a sphere of radius $R$ and, fixing its volume $V$, we increase its surface area $S$ to obtain a given reduced volume $v$, which is defined as the ratio between $V$ and the volume of a sphere with surface area $S$, $v=\frac{3\sqrt{4\pi}V}{S^{3/2}}$. For a sphere $v=1$ and $v<1$ otherwise. During the area increase, we solve the shape that minimizes the Helfrich energy. Once the prolate shape has been obtained, we initialize the lipid densities on each monolayer close to their equilibrium state for the given shape, i.e. satisfying $\rho^{\pm}=\rho_{0}(1\mp dH)$. To perturb the initial density profiles, we add a localized perturbation $\delta\rho^{\pm}=\delta\breve{\rho}^{\pm}(1\mp dH)$, where $\delta\breve{\rho}^{\pm}=\delta\breve{\rho}^{\pm}_{m}\,f(\theta,\phi)$ is the perturbation of the densities at the neutral surfaces of each monolayer, $\delta\breve{\rho}^{\pm}_{m}$ is the maximum value of the perturbation at the outer and inner monolayers respectively, and $f(\theta,\phi)$ is a function with values from 0 to 1 of the angles $(\theta,\phi)$ of a set of spherical coordinates adapted to the prolate shape. Following Dimova et al. (2006); Rahimi & Arroyo (2012), we choose $\kappa=10^{-19}$\text{,}\mathrm{J}$$, $k_{S}=5\times 10^{-2}$\text{,}\mathrm{J}\text{,}\mathrm{m}$^{-2}$, $b_{I}=10^{9}$\text{,}\mathrm{J}\text{,}\mathrm{s}\text{,}\mathrm{m}$^{-4}$, $\mu=5\times 10^{-10}$\text{,}\mathrm{J}\text{,}\mathrm{s}\text{,}\mathrm{m}$^{-2}$, and the dilatational viscosity $\lambda=0$ (this parameter seems to play a minor role in the dynamics).

Example 1: Relaxation dynamics of a density disturbance in an axisymmetric vesicle of 200 nm. To compare with (Rahimi & Arroyo, 2012), we start by examining a small vesicle ($R=200$ nm) with a reduced volume $v=0.99$, to which we apply a disturbance of $5\%$ in the outer monolayer, $\delta\breve{\rho}_{m}^{+}/\rho_{0}=5\%$, with a distribution $f(\theta)=\tanh\left((w-\theta)/\pi\right)$, where $\displaystyle w={\pi}/{10}$ controls the width of the disturbance. We show some snapshots of the dynamics along with the time-evolution of the dissipation and the main energy contributions, see figure 15. Again, we observe that the total energy $\mathcal{F}$ decays with time (figure 15A), as expected. Furthermore, from figure 15B we observe that the largest energetic component is $\mathcal{F}_{\text{H}}$, the Helfrich energy. However, it does not play a significant role in this problem since its variation is very small. Instead, we observe that the relaxation of the stretching energy in the upper monolayer, which transiently increases that of the lower monolayer, is the main driver of the dynamics (see figure 15B). In snapshot III, we can observe how the local density asymmetry results in a small but noticeable shape change, whose signature can be seen in the curvature energy. Note that, given the versatility of subdivision surfaces to deal with meshes of arbitrary connectivity, we have used a surface mesh with a much higher resolution at the pole where the density disturbance is imposed, see figure 15C.

These dynamics can be rationalized introducing several time-scales for this model following Rahimi & Arroyo (2012). Gradients of the average density relax with a time-scale given by $t_{4}=\mu/k_{S}$, as they are driven by stretching energy and dragged by shear dissipation. This time-scale is size-independent, and usually very fast, $t_{4}\approx$10\text{,}\mathrm{ns} for our choice of model parameters. Gradients of density differences between monolayers are also penalized by the stretching energy. However, at fixed shape, these gradients relax by intermonolayer slippage. Indeed, density differences have been shown to diffuse with a diffusivity $D=k_{S}/b_{I}$ (Evans & Yeung, [1994](#bib.bib38)), which results in a time-scale $t_{1}=\bar{S}/D=\bar{S}b_{I}/k_{S}$, where $\bar{S}$ is the area of the density disturbance. However, density differences can also relax by curving the membrane, not mobilizing intermonolayer slippage, with a time-scale given by $t_{2}=\sqrt{\bar{S}}\mu/(k_{S}d)$. For the $200$ nm vesicle, we find that $t_{1}\approx 0.15$1\text{\,}\mathrm{ms} and $t_{2}\approx$1\text{,}\mathrm{\SIUnitSymbolMicro s}. All these time-scales are apparent in figure [15](#S5.F15)A and highlight the dramatic gap between time-scales in this model, which need to be resolved by the simulations. To address this challenge, we adapt the time-step as shown in figure [15](#S5.F15)D, with time-steps spanning six orders of magnitude, from $0.1\text{\,}\mathrm{ns}$ to $0.1\text{\,}\mathrm{ms}$. To adapt the time-step we follow the following prescription: if Newton’s method is solved less than $N_{S}$ steps, with $N_{S}$ given initially (usually a number between 4 and 6), we increase $\Delta t^{n+1}=f\Delta t^{n}$ with $f$ a scaling factor greater than $1$. If, however, Newton’s method does not converge in $N_{S}$ steps, we reduce $\Delta t^{n+1}$ as $\Delta t^{n+1}=\Delta t^{n}/f$. This adaptive time-stepping algorithm allows us to perform the simulation in less than 300 time-steps, whereas a fixed time-step algorithm with the required initial resolution would need 10 million of time-steps. To show that the dynamics is not affected by the adaptive time-stepping, we plot the difference in the total energy between a simulation with a fixed and very small time-step ($\Delta t=$0.1\text{\,}\mathrm{ns}) and the simulation with the adaptive time-steping for the first 100 ns of dynamics, which shows a difference smaller than $0.1\%$ (figure 15E). Another important aspect of the numerical method is the global conservation of mass and volume. Conservation of the total mass depends on the local mass conservation imposed weakly through the process operator, whereas conservation of volume is imposed as a non-linear constraint at every time-step. We show the time-evolution of the relative error in total mass for the outer and inner monolayers in figure 15F, where we observe errors smaller than $10^{-2}\%$. We find errors in enclosed volume conservation smaller than $10^{-4}\%$.

Example 2: Relaxation dynamics of a density disturbance in a non-axisymmetric vesicle of 200 nm. To further show the versatility of the numerical method, we examine the dynamics of a non-axisymmetric system, in which the density disturbance is larger, $\delta\rho_{m}^{+}=25\%$, and not aligned with the symmetry axis of the prolate initial vesicle (see figure 16 and Movie 2). We observe a similar dynamics, now with a larger bulge due to the larger density difference, and with an initial stretching energy 4 times larger than the bending energy.

Example 3: Relaxation dynamics of a density disturbance in an initially axisymmetric vesicle of 2 micron. Finally, we analyze a vesicle of $R=$2\text{,}\mathrm{\SIUnitSymbolMicro m} with $\delta\breve{\rho}_{m}^{+}/\rho_{0}=5\%$. For this size, the stretching energy becomes even more dominant than for the $R=$200\text{\,}\mathrm{nm} vesicle. Indeed, the relative influence between the different energetic components is highly size-dependent. Given two vesicles, say 1 and 2, related by a geometric scaling factor $X$, we have that $\mathcal{F}_{\text{H}}(2)=\mathcal{F}_{\text{H}}(1)$ (the Helfrich energy is size independent), whereas $\mathcal{F}_{\text{S}}(2)=X^{2}\mathcal{F}_{\text{S}}(1)$. In agreement with this, the dynamics for $R=$2\text{,}\mathrm{\SIUnitSymbolMicro m} show the formation of a large bulge that affects the shape of the whole vesicle and with a stretching energy 20-fold larger than the Helfrich energy (see figure [17](#S5.F17)). The time-scales associated to this problem are $t_{1}\approx$15\text{\,}\mathrm{ms} and $t_{2}\approx$10\text{,}\mathrm{\SIUnitSymbolMicro s}, with $t_{4}=$20\text{\,}\mathrm{ns} as before. In agreement with these time-scales, we observe again the first energy decrease in a scale comparable with $t_{4}$, and a total duration of the relaxation dynamics of $10\text{\,}\mathrm{ms}$, similar to $t_{1}$. In figure 17B we plot the different dissipation contributions, shear viscosity and intermonolayer friction, in a log-log plot. This plot shows that, during the initial equilibration of the total density and during the bulge formation, shear dissipation dominates. However, at later stages, density differences relax due to intermonolayer slippage. In this time-adaptive simulation, the smallest and largest time-steps differ by 7 orders of magnitude.

Interestingly, in the initial stages of the bulge formation (figure 17III), we observe that a pattern resembling buckling forms at the edge of the bulge, presumably caused by a transient and local compression in a large enough region compared to the Föppl-von Kármán length-scale $\displaystyle l_{FvK}=\sqrt{\kappa/\sigma}\approx$5\text{,}\mathrm{nm}, where surface tension is dominated by stretching energy $\sigma=k_{S}\left(\left(\rho^{\pm}/\rho_{0}\right)^{2}-1\right)\approx 10^{-2}$ J$\bcdot$m*-2* for $\rho^{\pm}=1.05\rho_{0}$. This kind of transient buckling deformation is a three-dimensional phenomenon that could not develop in the axisymmetric simulations by Rahimi & Arroyo ([2012](#bib.bib90)). To examine this phenomenon further, we zoom figure [17](#S5.F17)III in the region where the pattern forms, see figure [18](#S5.F18). The formation of the pattern does not result from an increase of the total energy, which suggest that it is not caused by a numerical instability of our method. In Figs. [18](#S5.F18)I and II, we show the velocity field of the outer monolayer near the bulge at two different instants during the pattern formation. After the pattern has formed, Figs. [18](#S5.F18)III and IV, the amplitude of the bulge continues to increase, and the oscillatory deformation pattern progressively disappears, see figure [18](#S5.F18)V. The rest of the dynamics is similar to that obtained by Rahimi & Arroyo ([2012](#bib.bib90)), which suggest that the pattern forms due to an initial buckling instability that does not affect the final fate of the dynamics. To further test the stability of our scheme, we used a finer mesh and found the same dynamics; fluctuations develop with the same length scale, which suggest that it is a physical outcome of the model rather than an instability of the method. Our model lacks the dissipative forces induced by the bulk medium, which may modify this buckling-induced transient pattern formation. Indeed, the size of the disturbance is close to the Saffman-Delbrück length, $l_{\text{SD}}=$5\text{\,}\mathrm{\SIUnitSymbolMicro m}, and therefore bulk dissipation could start playing a role (Saffman & Delbrück, 1975; Arroyo & Desimone, 2009)

5.3 The cell cortex: A viscous layer driven by active tension

The elementary model of the actomyosin cortex introduced in section 3.3 exhibits a non-trivial phenomenology and reproduces to a large extent the mechanics of cells in different processes, such as during cytokinesis (Turlier et al., 2014) or in rheological assays (Torres-Sánchez, 2017). Here, we focus on the ability of this model to describe adhesion-independent cell migration.

In this kind of migration (Bergert et al., 2015; Ruprecht et al., 2015), cells develop a persistent cortical flow from the front to the rear of the cell that propel the cell forward by unspecific friction under confinement, see figure 19A. This friction is independent of specific adhesion molecules. Adhesion-independent locomotion plays a major role in three-dimensional cell migration through the extracellular matrix or in confined environments (Poincloux et al., 2011; Liu et al., 2015).

Adhesion-independent migration raises several questions. First, what is the mechanism by which cells acquire such a polarized state? Second, how can this flow be made persistent to allow for a self-sustained motion? And, how does the tight interplay between interfacial flows on the cortex and cell shape changes manifest itself in this process? Several models based on the theory of active gels have been developed over the past decade to try to answer these questions (Hawkins et al., 2011; Tjhung et al., 2012; Callan-Jones & Voituriez, 2013). In these models, myosin-mediated contraction of the cell cortex is identified as the main driver of the self-polarization. In particular, a spatial fluctuation in myosin activity can lead to a tension gradient in the cortex. This tension gradient triggers cortical flows, which further reinforce the gradient of tension, see figure 19B. This mechanism works against actin turnover, which tries to homogenize the system. Thus, adhesion-independent migration depends on the competition between myosin activity and actin turnover. Most of previous models rely on simple one-dimensional or fixed-shape assumptions that cannot address the effect of shape in locomotion. Recently, Callan-Jones et al. (2016) studied the shape transformations that cells suffer during migration, but this work is restricted to small deformations around a sphere. Here, we present, to our best knowledge, the first numerical results of a fully three-dimensional and nonlinear model connecting cortical flows and cell shape dynamics during locomotion. This work opens the door to a more systematic study of adhesion-independent migration in the future.

A critical ingredient controlling the formation of a self-polarized cortical flow is myosin activity, which is described by the function $\xi(\rho)$ in our model. If this function is constant, as assumed in previous works (Turlier et al., 2014), then the active tension is proportional to cortical thickness, $\gamma=\xi_{0}\rho$. In this case, the positive feedback illustrated in figure 19B leads to an instability with unbounded actin accumulation at the rear of the cell. Recently, Chugh et al. (2017) found that active tension does not depend linearly on cortical density in general. They found that in mitotic cells tension depends non-monotonically on cortical thickness, which they identified as a proxy for filament length. They proposed a conceptual model according to which active tension would be modulated by network architecture. Along the lines of this work, here we model the a dependence of specific contractility on cortical thickness as

[TABLE]

where $\xi_{0}$ measures a basal myosin activity and $\omega$ characterizes its dependence with cortex thickness. This leads to an active tension

[TABLE]

which has a maximum at $\rho=\omega\rho_{0}$; at equilibrium $\gamma_{0}=\gamma(\rho_{0})=\xi\rho_{0}(1-1/3\omega^{2})$. We note that the second term in the active tension looks very similar to the osmotic contribution introduced by Callan-Jones & Voituriez (2013) to stabilize the dynamics of polarization.

Following the experimental work by Ruprecht et al. (2015), we examine the migration of cells confined between two plates. To represent this confinement mathematically, we introduce a free energy contribution of the form

[TABLE]

where $z$ is a coordinate perpendicular to the plates, and $U(z)$ is a repulsive potential modelling contact with the plates and given by

[TABLE]

with $K_{c}$ and $\delta_{c}$ characterizing the strength and the width of the repulsive interaction respectively.

We now perform simulations of the model on a model cell of average radius $R=$5\text{,}\mathrm{\SIUnitSymbolMicro m}. Material parameters are obtained from literature $\rho_{0}=$500\text{\,}\mathrm{nm} (Clark et al., 2013), $\mu=10$\text{,}\mathrm{kPa}\text{,}\mathrm{s} (Bergert *et al.*, [2015](#bib.bib13)), $\tau=10$\text{\,}\mathrm{s} (Fritzsche et al., 2013), $\xi_{0}=1-10$\text{,}\mathrm{kPa} (Chugh *et al.*, [2017](#bib.bib25)). We also choose $\omega=2\sqrt{3}/3$. We first compress the cell between the plates with $h=$4\text{\,}\mathrm{\SIUnitSymbolMicro m} and let the system relax. To drive the cortex out of the equilibrium state at constant density $\rho=\rho_{0}$, we perturb the system with a gradient in density of $1\%$ in the $x$ direction simulating a possible fluctuation of myosin activity within the cortex, see figure 20A left. We first simulate the system with $\xi_{0}=$1\text{,}\mathrm{kPa}. For this small degree of contractility, the tension difference generated by the initial thickness perturbation is not high enough to overcome cortical turnover, and the system quickly relaxes to a situation of homogeneous cortical density, see figure [20](#S5.F20)A top right. For a higher value of myosin activity, $\xi_{0}=$10\text{\,}\mathrm{kPa}, the cortical flow generated by the activity gradient overcomes turnover, and the cell becomes self-polarized with a sustained cortical flow, see figure 20 bottom right. Together with the flow, the cell experiences a shape change during the transient dynamics towards the steady self-polarized state. Our ALE method is able to sustain such shape changes without remeshing, see figure 20B. More remarkably, we observe that the mesh is not affected by the constant flow of actin from the front to the rear of the cell. In a Lagrangian framework, such steady state would continuously distort the mesh, and very frequent remeshing would be required. Finally, we observe that, since tension is not a monotonic function of cortical thickness, it exhibits a maximum between the front and the rear of the cell, see figure 20C. This self-polarized state, however, cannot lead to cell migration by itself unless we introduce a mechanical interaction with the confining plates. To represent unspecific friction, we introduce the dissipation potential

[TABLE]

where $\eta_{c}$ measures friction with the plates and $U^{\prime}(z)$ identifies the pressure exerted by the cell on the plates. This pressure is equal to the internal pressure of the cell $P$, which is essentially determined by the cell radius of curvature and its surface tension and is $P\approx$0.3\text{,}\mathrm{kPa} in our simulations. Resorting to experimental measurements of the product of $\eta_{c}P=1-10^{4}$\text{\,}\mathrm{kPa}\text{\,}\mathrm{s}\text{\,}\mathrm{m}$^{-1}$ on somewhat larger cells (Bergert *et al.*, [2015](#bib.bib13)), we choose $\eta_{c}=$600\text{\,}\mathrm{s}\text{\,}\mathrm{m}$^{-1}$. We note, however, that our results are largely independent of friction because we do not consider a hydrodynamical resistive force in the relatively unconfined situation of cell motion between parallel plates. We repeat the previous simulation at $\xi_{0}=$10\text{\,}\mathrm{kPa}, and observe how the self-polarization of the cell now leads to cell migration, in a direction opposite to the cortical flow, see figure 21 and Movie 3. We note that, aside from a small disturbance of cortical velocity due to friction with the plates, the velocity field of actin in the steady state is the sum of a constant center of mass velocity plus a velocity profile similar to the one in figure 20A bottom right (data not shown). For these simulations and to deal with cell migration, we consider the following ALE parametrization

[TABLE]

where we impose zero net displacement due to the offset, $\int_{\Gamma_{t}}h(t)\bm{M}dS=\bm{0}$, and incorporate a rigid body translation $\bm{R}(t)$ as an unknown.

In conclusion, our theoretical and computational framework allows us to formulate and simulate thin and curved active gels with high generality. We have illustrated that this approach can be used to examine systematically adhesion-independent cell migration under confinement. Remarkably, our ALE formulation allows us to deal with the shape changes that the cell experiences during self-polarization and confinement, and with the steady cortical flows that are established.

6 Summary, discussion and future work

We have introduced a novel theoretical and computational framework to model and simulate fluid surfaces. Fluid surfaces are a common motif in cell and tissue biology. Thanks to increasingly quantitative biophysical experiments, there is a growing need for accurate theoretical predictions. Yet, modelling these systems requires overcoming significant theoretical and computational challenges, which we have addressed in this work. First, based on time-evolving parametrizations, we have rigorously extended the notion of ALE methods to fluid surfaces. We have also used Onsager’s formalism, a general variational framework for the dissipative dynamics of soft-matter systems, to derive thermodynamically consistent models of fluid surfaces coupling multiple physics in a fully geometrically non-linear manner. From a numerical perspective, we have proposed a new framework for the simulation of fluid surfaces based on a variational and nonlinearly stable time-integrator rooted in Onsager’s variational formalism, allowing us to bridge time-scales over 7 decades, and on a combination of subdivision and linear finite elements.

We have applied the previous theoretical and numerical methods to derive the governing equations and simulate the dynamics of canonical models of fluid surfaces with unprecedented generality (in three dimensions, for general shapes, and accounting for full geometric nonlinearity). We have first studied the dynamics of lipid bilayers in a number of interesting assays, including membrane relaxation, deflation due to osmotic shocks or perturbations due to density disturbances. Our framework opens new possibilities in the study of shape pattern formation under dynamical changes in lateral strain or osmotic conditions in supported membranes (Staykova et al., 2013) beyond axisymmetry, relevant to cell membrane mechano-adaptation (Kosmalska et al., 2015). Our method could also be useful to understand the effective rheology of a bilayer populated by transmembrane proteins, limiting inter-monolayer slippage in a heterogeneous manner, which could explain the unexpected and highly viscous behaviour of complex biomembranes (Campillo et al., 2013), or coupled to additional fields describing the concentration of membrane proteins to understand the dynamics of curvature sensing and generation (see (Arroyo et al., 2018; Baumgart et al., 2011) and references therein). While interfacial hydrodynamics are dominant at length-scales smaller than the Saffman-Delbrück length, the bulk hydrodynamics may be a relevant ingredient in processes involving larger scales. Including the bulk hydrodynamics is straightforward conceptually, but requires specialized computational methods, such as immersed boundary methods (Liu et al., 2006).

We have also applied our methodology to model and simulate the cell cortex. Our model is based on a viscous isotropic fluid layer, which is able to reproduce a number of rheological experiments and could be employed to infer material parameters in conjunction with experiments (Torres-Sánchez, 2017). Here, we have shown that our model is capable of reproducing adhesion-independent cell migration. Our simulations show how our ALE method can deal with the shape transformations that cells experience during migration and at the same time it can withstand steady flows from the front to the rear of the cell during migration. While our model for the cortex can reproduce a number of cellular behaviours, it is insufficient to reproduce phenomena where the transient elastic behaviour of the cortex becomes important, e.g. during laser ablation (Saha et al., 2016), or situations in which the orientational order of actin filaments becomes relevant (Reymann et al., 2016). This would require introducing tensorial fields on the surface (Nestler et al., 2018). Furthermore, a more detailed mechano-chemical model of activity, the explicit treatment of the cytosol, and models capable of spontaneously producing polarization would provide a more complete understanding of the mechanics of the cortex. These and other extensions of the active gel model presented here are enabled by the theoretical and computational tools introduced here.

Acknowledgements

We acknowledge the support of the European Research Council (CoG-681434), the European Commission (project H2020-FETPROACT-01-2016-731957), the Spanish Ministry of Economy and Competitiveness/FEDER (DPI2015-71789-R to MA and BES-2012-05489 to ATS), and the Generalitat de Catalunya (SGR-1471, ICREA Academia award to MA). We also thank Nikhil Walani, Sohan Kale and Daniel Santos-Oliván for useful discussions.

Appendix A Relation between Lagrangian and ALE velocities

[TABLE]

Appendix B Relation between Lagrangian and ALE time-derivatives

[TABLE]

Here we identify $\left[D\tilde{f}\tilde{\bm{c}}\right]\circ\bm{\psi}^{-1}$ as the pull-back of $\bnabla f\bcdot\bm{c}$, where $\bnabla f$ is the surface gradient of $f$.

Appendix C Rate-of-deformation tensor in terms of velocities

To obtain the form $\bm{d}$ in terms of $\bm{V}$, let us consider the components of $\bar{\bm{g}}$, which coincide with those of $\bm{g}$ in the convected basis by $\bm{\phi}$ given by the tangent vectors $\bm{e}_{a}=\partial_{a}\bm{\phi},~{}a=1,2$

[TABLE]

Then, we have

[TABLE]

where we have used the conmutativity of partial derivatives, the definition of covariant derivative

[TABLE]

the orthogonality of $\bm{N}$ to the tangent space of $\Gamma_{t}$ $\bm{e}_{a}\bcdot\bm{N}=0$, and the definition of the second fundamental form

[TABLE]

Appendix D Weak form of an inextensible monolayer with bending rigidity

To derive the weak form of the problem, we rewrite the material time derivative of the free energy (equation (61)) as

[TABLE]

where we have used that $\bnabla_{a}Hg^{ab}=\bnabla_{a}k^{ab}$, that

[TABLE]

and taken into account that the last two terms are null Lagrangians. Thus, variations of the velocity field around the solution $\bm{V}$ of the form $\bm{V}+\bm{U}$ lead to

[TABLE]

Equivalently, taking variations of the dissipation potential (equation (58)), we get

[TABLE]

where, again, we have set to zero null Lagrangians. Variations of the inextensibility constraint result in

[TABLE]

Finally, the last two terms have trivial variations

[TABLE]

and

[TABLE]

Collecting all these variations, we have the following statement of stationarity

[TABLE]

which should hold for all admissible variations $\bm{U}$, where

[TABLE]

from where one can identify

[TABLE]

with

[TABLE]

and

[TABLE]

Finally,

[TABLE]

Appendix E Weak form of the three-dimensional non-linear Seifert-Langer model

In this case, we focus on the stretching energy, dilatation dissipation and intermonolayer friction, since the rest of terms were already derived for an inextensible monolayer (see D). The rate of change of the stretching energy is

[TABLE]

Accounting for multiplicative factors, term 1 leads to

[TABLE]

and term 5 to

[TABLE]

Summing them, we note that their first terms cancel out with each other. Rearranging the last terms, we get

[TABLE]

plus null Lagrangians, which we neglect for the sake of simplicity since we are dealing with a closed surface. Let us define the stress tensors

[TABLE]

[TABLE]

[TABLE]

and the normal stress vector

[TABLE]

Then, equation (167) can be rewritten as

[TABLE]

Terms 2 plus 4 lead to

[TABLE]

Term 3, neglecting null Lagrangians,

[TABLE]

Altogether, we can write the rate of change of the free energy as

[TABLE]

Thus,

[TABLE]

and

[TABLE]

From variations of the dilatation dissipation potential, we get

[TABLE]

and

[TABLE]

Finally, variations of the intermonolayer friction dissipation potential lead to

[TABLE]

Appendix F Discrete free energy and dissipation potentials for an inextensible viscous monolayer with bending energy

We have defined the discrete Helfrich energy,

[TABLE]

where we have split integration on $\Gamma_{t}$ as a sum of integration on the curved triangles $\Gamma_{t}^{E}$, which are evaluated at the parametric domains $\tilde{\Gamma}$. Functions $H\left(\mathsf{h}\right)$ and $J\left(\mathsf{h}\right)$ describe the mean curvature and the surface Jacobian in terms of the discretized parametrization; these can be computed by plugging the form of $\bm{\psi}$ (equation (122)) in terms of $\mathsf{h}$ in the expressions for the curvature and metric in the natural or convected basis of the parametrization. We have also defined the matrices representing dissipation and the inextensibility constraint

[TABLE]

where, in the last two equations, the functions $N_{I}^{E}$, interpolating the surface tension $\gamma$ are composed with $O^{-1}$ and evaluated at the macroelement. We note that $\bnabla$ is the covariant derivative, calculated from partial derivatives in parametric space and using Christoffel symbols (Do Carmo, 1992). We also note that we have also discretized the rate of change of volume$\dot{\Omega}$ as $(\Omega-\Omega^{n})/\Delta t^{n}$ instead of discretizing equation (65) directly, similarly to our discretization of the energy release rate. This leads to a discrete dynamics that keeps a constant volume by construction, up to numerical error, regardless of the value of $\Delta t^{n}$. To exercise this formulation, we compute $V(\mathsf{h})$ using Gauss theorem on the surface

[TABLE]

We finally note that we use Gauss quadrature in the reference element $\tilde{\Gamma}$, although other integration schemes specially suited for subdivision surfaces have been recently proposed (Jüttler et al., 2016).

Appendix G Discretization of mass conservation

We consider an implicit Euler scheme to discretize in time the process operator in the transport problem as in equation (115), which leads to

[TABLE]

In this case, $\bm{d}$ and $\bm{c}$ depend on $(\mathsf{h},\mathsf{a},\mathsf{b})$, but we do not write it for simplicity. This is a reaction-advection problem in $\rho^{n+1}$ and its discretization with finite elements has to be carefully considered, since Garlerkin methods cannot deal with large convective terms. Discretizing we obtain

[TABLE]

To deal with the convective term appropriately, we use the test functions

[TABLE]

following a Petrov-Garlerkin method in which the weight functions do not coincide with the basis functions used in the approximation of the solution $\rho{n+1}$. This method is called stream-upwind Petrov Garlerkin (SUPG) (Donea & Huerta, 2003), which is able to treat the convective term of the transport problem by adding numerical diffusion controlled by the SUPG parameter $\gamma_{s}$. Then, the weak form is

[TABLE]

where here $J$ is also a function of $\mathsf{h}$. This equation can also be written as a linear system

[TABLE]

with

[TABLE]

and

[TABLE]

We note that $\hat{\mathsf{M}}$ is not symmetric and $\hat{\mathsf{M}}$ and $\hat{\mathsf{s}}$ depend non-linearly on $(\mathsf{h},\mathsf{a},\mathsf{b})$ through $\bm{d}$, $\bm{c}$ and $J$. The coupled system of finite element equations involving balance of linear momentum and mass transport, corresponding to the spatial discretization of equation (117), are solved simultaneously using a Newton-Raphson method.