Computational modeling of active deformable membranes embedded in 3D flows

TL;DR

This paper introduces a computational algorithm that models the dynamic behavior of active, deformable membranes in 3D flows by integrating active gel theory with elastic shell mechanics and fluid coupling techniques.

Contribution

The authors develop a novel computational method combining active gel theory, thin shell mechanics, and fluid dynamics to simulate active membranes of arbitrary shape in 3D flows.

Findings

Validated the algorithm with known Green's functions for cylindrical membranes.

Predicted a new non-axisymmetric instability in active membranes.

Demonstrated the method's application to cell deformation in shear flow.

Abstract

Active gel theory has recently been very successful in describing biologically active materials such as actin filaments or moving bacteria in temporally fixed and simple geometries such as cubes or spheres. Here we develop a computational algorithm to compute the dynamic evolution of an arbitrarily shaped, deformable thin membrane of active material embedded in a 3D flowing liquid. For this, our algorithm combines active gel theory with the classical theory of thin elastic shells. To compute the actual forces resulting from active stresses, we apply a parabolic fitting procedure to the triangulated membrane surface. Active forces are then dynamically coupled via an Immersed-Boundary method to the surrounding fluid whose dynamics can be solved by any standard, e.g. Lattice-Boltzmann, flow solver. We validate our algorithm using the Green's functions of [Berthoumieux et al. New J. Phys. 16, 065005 (2014)] for an active cylindrical membrane subjected (i) to a locally increased active stress and (ii) to a homogeneous active stress. For the latter scenario, we predict in addition a so far unobserved non-axisymmetric instability. We highlight the versatility of our method by analyzing the flow field inside an actively deforming cell embedded in external shear flow. Further applications may be cytoplasmic streaming or active membranes in blood flows.

Figures

Equations

$$\frac{{\,\mathrm{d}{f}}}{{\,\mathrm{d}{t}}}=\frac{\partial f}{\partial t}+\bm{v}\cdot\mbox{\boldmath$\nabla$}_{\bm{r}}f+\bm{F}\cdot\mbox{\boldmath$\nabla$}_{\bm{p}}f=\Omega$$

$$f_i({\mathbf{x}}_j+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i({\mathbf{x}}_j,t)+{\Omega }_i({\mathbf{x}}_j,t)\Delta t.$$

$$\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu \Delta \mathbf{u}+\frac{1}{\rho }\int \mathbf{f}\left({\mathbf{X}}^{\prime },t\right)\delta \left({\mathbf{X}}^{\prime }-\mathbf{x}\right)\mathrm{d}{}^2{X}^{\prime },$$

$${\mathbf{F}}_j$$

$$f_i({\mathbf{x}}_j+{\mathbf{c}}_i\Delta t$$

$$\frac{\mathrm{d}{\mathbf{x}}_n}{\mathrm{d}t}=\mathbf{u}\left({\mathbf{x}}_n(t),t\right)=\int \mathbf{u}(\mathbf{x},t)\delta \left(\mathbf{x}-{\mathbf{x}}_n\right)\mathrm{d}{}^{3}x.$$

$$w^{\text{SK}}(I_1,I_2)=\frac{{\kappa }_{\text{S}}}{12}\left((I_1^2+2I_1-2I_2)+C{I}_2^2\right).$$

$$w^{\text{HF}}=2{\kappa }_{\text{B}}(H-H_0{)}^2+{\kappa }_{K}K$$

$$-f_j=\left({\sigma }_{ij}^{out}-{\sigma }_{ij}^{in}\right)n_i$$

$$-\mathbf{f}=-f^{\alpha }{\mathbf{e}}_{\alpha }-f^n\mathbf{n}.$$

$${\nabla }_{\alpha }^{\prime }{t}_{\text{a}}^{\alpha \beta }+C_{\alpha }^{\prime \beta }{t}_{n\text{a}}^{\alpha }+f_{\text{e}}^{\beta }$$

$${\nabla }_{\alpha }^{\prime }{t}_{n\text{a}}^{\alpha }-C_{\alpha \beta }^{\prime }{t}_{\text{a}}^{\alpha \beta }+f_{\text{e}}^n$$

$${\nabla }_{\alpha }^{\prime }{m}_{\text{a}}^{\alpha \beta }+C_{\alpha }^{\prime \beta }{m}_{n\text{a}}^{\alpha }$$

$${\nabla }_{\alpha }^{\prime }{m}_{n\text{a}}^{\alpha }-C_{\alpha \beta }^{\prime }{m}_{\text{a}}^{\alpha \beta }$$

$$f_{\text{a}}^{\beta }$$

$$f_{\text{a}}^n$$

$${\mathbf{e}}_{\xi }=\frac{{\mathbf{x}}_1-\left({\mathbf{x}}_1\cdot \mathbf{n}\right)\mathbf{n}}{|{\mathbf{x}}_1-\left({\mathbf{x}}_1\cdot \mathbf{n}\right)\mathbf{n}|}.$$

$${\mathbf{e}}_{\eta }=\frac{\mathbf{n}\times {\mathbf{e}}_{\xi }}{|\mathbf{n}\times {\mathbf{e}}_{\xi }|}.$$

$${\bar{r}}_i(\xi ,\eta )\approx r_i^c+A_i\xi +B_i\eta +\frac{1}{2}\left(C_i{\xi }^2+D_i{\eta }^2+2E_i\xi \eta \right),$$

$${\chi }_i^2=\sum _{\nu =1}^N{\left(r_i^{\nu }-{\bar{r}}_i^{\nu }\right)}^2,$$

$$A_i={\bar{r}}_{i,\xi }\text{leavevmode}\text{leavevmode}\text{leavevmode}{B}_i={\bar{r}}_{i,\eta }\text{leavevmode}\text{leavevmode}\text{leavevmode}{C}_i={\bar{r}}_{i,\xi \xi }\text{leavevmode}\text{leavevmode}\text{leavevmode}{D}_i={\bar{r}}_{i,\eta \eta }\text{leavevmode}\text{leavevmode}\text{leavevmode}{E}_i={\bar{r}}_{i,\xi \eta }$$

$$\mathbf{A}={\bar{\mathbf{r}}}_{,\xi }{|}_{{\mathbf{r}}_c}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathbf{B}={\bar{\mathbf{r}}}_{,\eta }{|}_{{\mathbf{r}}_c}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathbf{C}={\bar{\mathbf{r}}}_{,\xi \xi }{|}_{{\mathbf{r}}_c}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathbf{D}={\bar{\mathbf{r}}}_{,\eta \eta }{|}_{{\mathbf{r}}_c}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathbf{E}={\bar{\mathbf{r}}}_{,\xi \eta }{|}_{{\mathbf{r}}_c}.$$

$$g_{\tilde{\alpha }\tilde{\beta }}=\left(\begin{array}{cc}\mathbf{A}\cdot \mathbf{A}& \mathbf{A}\cdot \mathbf{B}\\ \mathbf{B}\cdot \mathbf{A}& \mathbf{B}\cdot \mathbf{B}\end{array}\right),$$

$$C_{\tilde{\alpha }\tilde{\beta }}=\left(\begin{array}{cc}-\mathbf{C}\cdot \mathbf{n}& -\mathbf{E}\cdot \mathbf{n}\\ -\mathbf{E}\cdot \mathbf{n}& -\mathbf{D}\cdot \mathbf{n}.\end{array}\right),$$

$$g_{\xi \xi ,\xi }={\left(\mathbf{A}\cdot \mathbf{A}\right)}_{,\xi }=2\mathbf{A}\cdot \mathbf{C},$$

$${\bar{t}}_{\text{a}}^{\alpha \beta }(\xi ,\eta )=t_{\text{a}}^{c\alpha \beta }+A_{\text{a}}^{\alpha \beta }\xi +B_{\text{a}}^{\alpha \beta }\eta +\frac{1}{2}\left(C_{\text{a}}^{\alpha \beta }{\xi }^2+D_{\text{a}}^{\alpha \beta }{\eta }^2+2E_{\text{a}}^{\alpha \beta }\xi \eta \right)$$

$${\chi }_{\text{a}}^2=\sum _{\nu =1}^N(t_{\text{a}}^{\tilde{\alpha }\tilde{\beta }\nu }-{\bar{t}}_{\text{a}}^{\tilde{\alpha }\tilde{\beta }\nu }{)}^2$$

$$f_{\text{a}}^{\tilde{\beta }}$$

$$f_{\text{a}}^n$$

$$A_{\text{Voronoi}}=\frac{1}{8}\sum _{\nu =1}^N\left(\cot ({\alpha }_{\nu })+\cot ({\beta }_{\nu })\right)|{\mathbf{x}}_{\nu }|$$

$$t_{\text{a}\alpha }^{\text{leavevmode}\text{leavevmode}\beta }=\left(\begin{array}{cc}{t}_{\text{a}1}^{\text{leavevmode}\text{leavevmode}1}& t_{\text{a}1}^{\text{leavevmode}\text{leavevmode}2}\\ t_{\text{a}2}^{\text{leavevmode}\text{leavevmode}1}& t_{\text{a}2}^{\text{leavevmode}\text{leavevmode}2}\end{array}\right)$$

$$t_{\text{a}}^{\tilde{\alpha }\tilde{\beta }}=\mathfrak{D}\text{leavevmode}{t}_{\text{a}}^{\alpha \beta }\text{leavevmode}{\mathfrak{D}}^{-1},$$

$$t_{\text{a}\alpha }^{\text{leavevmode}\beta }=\left(\begin{array}{cc}{T}_{\text{a}}+T_{\text{a}}^z\delta (z)& 0\\ 0& T_{\text{a}}+T_{\text{a}}^{\varphi }\delta (z)\end{array}\right)$$

Full text

Computational modeling of active deformable membranes embedded in three-dimensional flows

Christian Bächer christian.baecher@uni-bayreuth.de Biofluid Simulation and Modeling, Theoretische Physik VI, Universität Bayreuth, Bayreuth, Germany

Stephan Gekle stephan.gekle@uni-bayreuth.de Biofluid Simulation and Modeling, Theoretische Physik VI, Universität Bayreuth, Bayreuth, Germany

Abstract

Active gel theory has recently been very successful in describing biologically active materials such as actin filaments or moving bacteria in temporally fixed and simple geometries such as cubes or spheres. Here we develop a computational algorithm to compute the dynamic evolution of an arbitrarily shaped, deformable thin membrane of active material embedded in a 3D flowing liquid. For this, our algorithm combines active gel theory with the classical theory of thin elastic shells. To compute the actual forces resulting from active stresses, we apply a parabolic fitting procedure to the triangulated membrane surface. Active forces are then dynamically coupled via an Immersed-Boundary method to the surrounding fluid whose dynamics can be solved by any standard, e.g. Lattice-Boltzmann, flow solver. We validate our algorithm using the Green’s functions of [Berthoumieux et al. New J. Phys. 16, 065005 (2014)] for an active cylindrical membrane subjected (i) to a locally increased active stress and (ii) to a homogeneous active stress. For the latter scenario, we predict in addition a so far unobserved non-axisymmetric instability. We highlight the versatility of our method by analyzing the flow field inside an actively deforming cell embedded in external shear flow. Further applications may be cytoplasmic streaming or active membranes in blood flows. Published as: Phys. Rev. E 99(6), 062418 (2019)

1 Introduction

Many biological processes such as cell division or locomotion depend on the ability of living cells to convert chemical energy into mechanical work [1]. A prominent mechanism to achieve such a conversion are motor proteins which perform work through a relative movement of cross-linked cytoskeletal filaments. This movement induces active stresses in the cell cortex [2, 3, 4, 5, 6, 7] which are transmitted via anchor proteins to the plasma membrane separating the interior of the cell from its surroundings. Active stresses have an inherent non-equilibrium character and are the basis of physically unique processes in active fluid layers such as instabilities [8, 9], the emergence of spontaneous flows [10] or the creation of geometrical structures [11, 12]. They are furthermore responsible for large shape deformations during cell morphogenesis [13, 14, 15], cell division [16, 17, 18, 19, 20, 21], cell locomotion [22, 23, 24, 25, 26], cell rheology [27, 28] or spike formation on artificial vesicles [29].

In recent years a theoretical framework has been developed describing cytoskeleton and motor proteins together as an active continuous gel [30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. This active gel theory in general treats the cell cortex as a viscoelastic medium with additional active contributions. On time scales short compared to the viscoelastic relaxation time active gel theory can be formulated in the elastic limit [3]. For thin active 2D membranes embedded in 3D space, active gel theory can be reformulated into force balance equations using the formalism of differential geometry [3, 6]. Any active stress in the membrane is then balanced by a counter-stress, usually due to viscous friction, from the external medium. The force balance equations can be applied on fixed membrane geometries such as spheres, cylinders or flat layers. This often results in the prediction of regions in parameter space where the prescribed shape is expected to become unstable [16, 40, 29, 11, 41, 42]. Despite being a powerful qualitative tool, such calculations cannot make statements about the precise shape of the active membrane after the instability has set in.

To obtain actual shape predictions, a number of works start instead from a parametrized, free membrane shape which is adjusted so as to fulfill the force balance equations for a given set of parameters and boundary conditions. This procedure can be carried out either analytically [43, 3] or numerically [20, 21, 44, 15, 45, 46]. For example, Berthoumieux et al. [3] derived Green’s functions to predict the deformation of an infinitely long cylindrical active, elastic membrane resulting from the application of a point active stress. For certain parameter ranges, these Green’s functions exhibit divergences which have been interpreted as mechanical (buckling or Rayleigh-Plateau-like) instabilities of the cylindrical membrane. Sain et al. [21] investigated the axisymmetric dynamics of the furrow in cytokinesis. Turlier et al. [20] computed the time evolution of an axisymmetric membrane undergoing cytokinesis by advecting tracer points discretizing the membrane. Callan-Jones et al. [15] predicted a transition to a polarized cell shape because of an instability of the cell cortex. Reymann et al. [45] matched axisymmetric, theoretical results to observed cell shapes during cytokinesis using measured velocity and order parameter fields as input for the theory. Heer et al. [47] determined the equilibrium shape of an elastic tissue layer folded into a deformable ellipsoidal shell, where myosin activity is incorporated as a preferred curvature of the shell. Klughammer et al. [48] analytically calculated the flow inside a sphere that slightly deforms due to a traveling band of surface tension, which mimics cortical active tension, under the assumption of rotational symmetry. In all of these works the full dynamics of the external fluid was neglected and nearly all of them restricted themselves to deformations from simple rest shapes in the steady state. There currently exists no analytical or numerical method to compute the dynamical deformation of an arbitrarily shaped active membrane immersed in a 3D moving Newtonian fluid.

In this work, we develop a computational algorithm to predict the dynamic deformation of arbitrarily shaped thin membranes discretized by a set of nodes connected via flat triangles. Starting from prescribed active stresses on the discretized membrane our algorithm computes the corresponding forces on every membrane node via a parabolic fitting procedure taking into account the full deformed surface geometry. Knowledge of the nodal forces then enables the dynamic coupling to a surrounding 3D fluid via the Immersed-Boundary Method (IBM). The Navier-Stokes dynamics for the surrounding fluid is solved here by the Lattice-Boltzmann method (LBM), but other flow solvers can straightforwardly be incorporated. With this, our method allows to study the dynamic evolution of active membranes in general external flows. It thus builds a bridge between the extensive literature on active fluids and the similarly extensive work on elastic cells, vesicles and capsules in flows [49, 50, 51, 52, 53, 54, 55, 56]. In biological situations often a flowing environment is present rendering the dynamic coupling of external fluid and membrane deformations necessary. Our proposed method allows such a coupling and thus the computation of dynamically evolving non-equilibrium shapes. Possible applications of our method include the study of active cell membranes inside the blood stream or active cellular compartments in cytoplasmic streaming flows.

First in section 2, we extensively describe the LBM and IBM for a dynamic coupling of an elastic membrane and a suspending fluid. In section 3 we start with the problem formulation in the framework of thin shell theory using differential geometry. Next, we describe our algorithm for three dimensional active force calculation in section 4. In section 5 we provide an in-depth validation of our algorithm based on the analytical results for a cylindrical active membrane in case of infinitesimal deformations by Berthoumieux et al. [3]. As an application, section 6.1 presents the flow field inside a dividing elastic cell, where active stresses trigger membrane deformations that in turn lead to fluid flow. In section 6.2 we analyze the same system in externally driven shear flow. Eventually, we conclude our work in section 7.

2 Dynamic coupling of membrane and fluid

2.1 Lattice-Boltzmann method for fluid dynamics

The Lattice-Boltzmann method (LBM) is an efficient and accurate flow solver which is well described in the literature [57, 58, 59, 60]. In the following we therefore summarize only the basic concepts. In contrast to macroscopic, e.g. finite element, methods based directly on the discretized Navier-Stokes equation (NSE), LBM starts from the Boltzmann equation which is a common tool in statistical mechanics. Using Chapman-Enskog analysis the NSE are recovered from the Boltzmann equation [61, 60].

The Boltzmann equation provides insight into a system on mesoscopic length scales by means of the continuous particle distribution function $f(\bm{r},\bm{p},t)$ where $\bm{r}$, $\bm{p}$ and $t$ refer to position of the particles, momentum of the particles and time, respectively. The expression $f(\bm{r},\bm{p},t){\,\mathrm{d}{\bm{r}}}{\,\mathrm{d}{\bm{p}}}{\,\mathrm{d}{t}}$ gives the probability to find a particle (fluid molecule) in the phase-space volume ${\,\mathrm{d}{\bm{r}}}{\,\mathrm{d}{\bm{p}}}$ at $(\bm{r},\bm{p})$ in the time interval $t$ to $t+dt$. The dynamic evolution of the particle distribution function $f(\bm{r},\bm{p},t)$ is given by

[TABLE]

with $\Omega$ being the collision operator accounting for re-distribution of molecules due to collisions.

Discretization of space, momentum space and time leads from the Boltzmann equation (1) to the Lattice-Boltzmann equation. The discretization of the spatial domain is carried out using a cubic Eulerian grid. The distance between the fluid nodes is $\Delta x=1$ in simulation units. In contrast to other methods, LBM also discretizes the momentum (velocity) space, i.e., $f=f(\bm{x}_{j},\bm{p}_{i},t)$ such that only a discrete set of velocities is allowed at each node. Here, we use the common D3Q19 scheme with 19 discrete velocity vectors, which is illustrated in figure 1 a). Thus, each node contains one population $f(\bm{x}_{j},\bm{p}_{i},t)$ for each momentum $\bm{p}_{i}$, which moves with the corresponding velocity $\bm{c}_{i}$ away from the node $\bm{x}_{j}$ within one time step. As an abbreviation the distribution functions are labeled by an index according to their discrete momentum/velocity, i.e., $f(\bm{x}_{j},\bm{p}_{i},t)=f_{i}(\bm{x}_{j},t)$. Finally, discrete time steps from $t$ to $t+\Delta t$ are considered. Under discretization the Boltzmann equation in (1) becomes the Lattice-Boltzmann equation [60]

[TABLE]

The numerical integration of the Lattice-Boltzmann equation in time is split into two steps, collision and propagation. Collision is done by an approximation of the collision operator in equation (2). The idea for the approximation of the collision operator is that the populations should relax towards the Maxwellian equilibrium distribution function in absence of driving forces. Here, we use the multiple relaxation time scheme. In the framework of the multiple relaxation time scheme the relaxation is done in moment space of the populations with individual relaxation rates for the different moments [58, 60].

Moments corresponding to conserved density and momentum do not have a relaxation time, while two relaxation rates related to bulk and shear viscosity, respectively, are chosen for moments corresponding to the stress tensor. For further discussion on the relaxation rates and moments used in our LBM implementation we refer to refs. [62, 58].

In the second step, the streaming, the populations after collision propagate according to the associated velocities.

We note that our LBM implementation supports adding thermal fluctuations to the fluid dynamics mimicking a given thermal energy $k_{B}T$ corresponding to fluctuating hydrodynamics [62]. Thermal fluctuations are taken into account by adding a random noise to those relaxed moments of the multiple relaxation time scheme which correspond to elements of the stress tensor [62, 58]. We note that our approach using thermal fluctuations is different from considering a separate temperature field [63, 64, 65] as e.g. required for convection flows [63]. Although they are not per se necessary, we employ thermal fluctuations in the present manuscript to speed up the onset of instabilities.

Although we do not consider solid walls in the present manuscript, we note that they can be realized by bounce back boundary conditions [60], where populations streaming towards nodes beyond a boundary are simply reflected. In case of moving elastic objects (2.2) having close contact to solid walls as it is considered e.g. in ref. [66], at least one fluid node between the object and the solid boundary is required in our method.

Typical LBM grids used in this work have dimensions of 126x72x72 or 100x100x100. Except in the case of the shear flow with moving boundaries in section 6.2, our LBM simulations use periodic boundary conditions for the fluid. A typical simulation run, e.g. in figure 6 consists of $10^{7}$ time steps. Onset of Rayleigh-Plateau like instabilities in figure 7 typically appears after approximately $8\times 10^{6}$ time steps. Simulations are performed with the simulation package ESPResSo [67, 68, 69] which has been extended to include thin membranes using the Immersed Boundary Method described next.

2.2 Immersed-Boundary method for membrane dynamics

The framework of the Immersed-Boundary Method (IBM) [70, 71, 72] allows for the coupling of a cellular membrane to the suspending fluid, where the fluid is simulated using LBM as described in the previous section 2.1. The IBM consists of two central steps: elastic forces acting on the membrane are transmitted to the fluid and due to no-slip boundary condition the mass-less membrane is advected with the local fluid velocity.

The membrane of a cell is represented by an infinitely thin elastic sheet in the framework of thin shell theory [73, 74, 75]. In our numerical simulations the membranes is discretized by a set of nodes that are connected by flat triangles [76, 51, 77, 78]. These represent the elastic membrane as a Lagrangian mesh immersed into the Eulerian fluid mesh as illustrated in figure 1 b).

Physical behavior of the membrane is characterized by appropriate constitutive laws which are described in section 3. The resulting forces are calculated on the deformed Lagrangian membrane mesh as described in section 4. To transmit these forces into the fluid, the incompressible NSE for the fluid velocity $\bm{u}(\bm{x},t)$ becomes [77]

[TABLE]

with $p$ the pressure, $\nu$ the kinematic viscosity, $\rho$ the density, and $\bm{f}$ the membrane force per area acting on the membrane, which is parametrized by $\bm{X}^{\prime}$. The Dirac $\delta$ function ensures that the forces only act at the position of the actual membrane. Correspondingly, the Lattice-Boltzmann equation (2) changes to

[TABLE]

After an update of the fluid dynamics (LBM algorithm), which now include the forces from the membrane $\bm{F}_{j}$, the mass-less membrane nodes are advected with the local fluid velocity thus satisfying exactly the no-slip boundary condition [71]. Moving with the local fluid velocity is expressed for a membrane node $\bm{x}_{n}\in\bm{X}^{\prime}$ by [70, 71]

[TABLE]

with $\bm{u}\left(\bm{x}_{n},t\right)$ being the fluid velocity at the position of the membrane node. In simulations equation (6) is integrated numerically using an Euler scheme in order to move each membrane node from time $t$ to $t+\Delta t$.

The core of both steps, force transmission and movement with local fluid velocity, is the interpolation between the Eulerian fluid grid and the Lagrangian membrane grid. Considering the discretization and the resulting spatial mismatch of membrane nodes and fluid nodes, as illustrated in figure 1 b), it becomes clear that an interpolation is necessary. On the one hand the ideal point force at the site of a membrane node must be spread to the adjacent fluid nodes. On the other hand to obtain the local fluid velocity at the site of a membrane node the velocity of the adjacent fluid nodes is interpolated. This means the $\delta$ distributions in equation (4) and (6) must be discretized. The interpolation between fluid and membrane is carried out by an eight-point stencil. As illustrated for two membrane nodes by the blue shaded cubes in figure 1 b), for each membrane node a cube containing the eight nearest fluid nodes is considered. A linear interpolation between the eight fluid points is performed [79].

A typical mesh for a cylinder in the present study contains 17100 Lagrangian membrane nodes and 34200 triangles. Simulating a cylindrical membrane, as it is partly shown in figure 4, is done using periodic boundary conditions for the fluid in all directions and for the membrane in the direction of the cylinder axis. The latter is achieved by connecting membrane nodes at the end of the box to those at the beginning of the box.

2.3 Validation of LBM/IBM for passive, elastic membranes

Our implementation of the force calculation and the LBM/IBM have extensively been tested for passive, elastic membranes in previous publications [52, 78, 80, 56, 81, 82].

The algorithm for the passive, elastic force calculation has been validated in references [52, 78] by comparison with exact results in static situations as well as for a capsule in shear flow. For red blood cells flowing through a rectangular channel very good agreement with in vitro experiments has been found [80, 56].

The mixed LBM/IBM has been extensively validated for suspensions of red blood cells and rigid particles in complex geometries. In [81] the concentration profile of cells across the channel diameter (affected by cross streamline migration, which is triggered by the passive elasticity of the cells) has been successfully compared to well-established literature results. In ref. [82] we have performed a validation based on the Zweifach-Fung effect for red blood cell suspension in branching channels. Furthermore, in the Supplementary Information of ref. [82] we have shown red blood cell shapes for a single cell in tube flow and in a rectangular channel, respectively, that are in very good agreement with previous studies using dissipative particle dynamics [83] and boundary integral method [56], respectively. The stability and behavior of stiff particles realized by IBM has been shown and validated in ref. [81] and [82] for a spherical particle based on the Stokes relation (sphere pulled through a quiescent fluid as well as a fixed sphere in homogeneous flow) and an ellipsoidal particle rotating in shear flow.

In addition to the well-established passive elasticity, in this paper we introduce active elastic forces into simulations of membranes, see section 4. We refer to section 5 for the validation in case of an active, elastic cell membrane.

3 Force balance in thin shell formulation including active stresses

3.1 Physical model

We consider the plasma membrane and the cytoskeletal filaments within the cell cortex, sketched in figure 2, as a single physical entity which for simplicity we denote in the following as ”the membrane”. Since this membrane is small in height compared to the cell diameter it is described as an infinitely thin shell. Thin shell theory treats the membrane as a two dimensional manifold embedded in the three dimensional environment, e.g., the intra- and extra-cellular fluid, thereby accounting for membrane curvature. For a condensed introduction into the required differential geometry on such manifolds as well as present conventions we refer the reader to ref. [84] and appendix A1 of this work, respectively. Here, we only note that Greek indices refer to coordinates on the membrane, i.e., $\alpha,\beta=1,2$, and Einstein summation convention is used. In this work, we consider the short-time limit of a purely elastic membrane noting that our algorithm can, without substantial difficulty, be extended to viscous or viscoelastic membranes.

The key quantity in the framework of thin shell theory are the in-plane surface stresses $t_{\alpha\beta}$ and moments $m_{\alpha\beta}$ [73, 74, 85]. The in-plane surface stress tensor or stress resultant [74] $t_{\alpha\beta}$ is the 3D elastic stress tensor for the material forming the membrane projected onto the membrane and integrated over the membrane thickness $h$ [74] (dimensions of a force per length, i.e., N/m) 111The term in-plane surface stresses mainly follows the terminology of Deserno [84]. Guven [86] and Deserno [84] denote the corresponding force vector with its in-plane components being $t_{\alpha\beta}$ as surface stress. For the same quantity Berthoumieux et al. [3] and Salbreux and Jülicher [6] as well as Guckenberger and Gekle [78] use the term in-plane tension tensor. Daddi-Moussa-Ider et al. [87] call this quantity stress tensor. Finally, Green and Adkins [74] call it stress resultant.. It is a function of in-plane coordinates $\alpha$, $\beta$ and contains only in-plane components such that its matrix representation has dimensions 2x2. Besides the in-plane surface stresses, we introduce the normal surface stress which is sometimes denoted a shearing force [73] or transverse shear surface stress $t_{n}^{\alpha}$ [76, 78]. The in-place surface stresses $t_{\alpha\beta}$ contain a contribution from passive, elastic stresses as well as from active, force-generating mechanisms. Both contributions superpose linearly and thus can be treated separately [3, 6]. Passive elastic stresses can be further split into different contributions such as, e.g., shear and bending resistances [50, 88, 78]. The moment tensor or stress couples [74] $m_{\alpha\beta}$ accounts for stress distribution across the membrane [74] (dimensions of a torque per unit length, i.e., N). We denote the corresponding normal surface moment by $m_{n}$. For explicit materials building up the membrane corresponding constitutive laws are required to derive explicit forms for in-plane surface stresses and moments.

Considering the passive, elastic properties of membranes in most cases constitutive laws are formulated in terms of a strain energy functional. For many biological membranes such as red blood cells the passive, elastic forces arise from the resistance to shear deformation of the cortical, cytoskeletal network and the resistance to bending deformation and area dilatation of the plasma membrane. For shear resistance and area dilatation Skalak et al. [89] proposed an appropriate energy density functional

[TABLE]

The energy density depends on the invariants $I_{1},I_{2}$ of the transformation between undeformed and deformed membrane [89, 50]. Both invariants are defined in equation (A2.3) and (A2.4) of the appendix. Resistance to shear is characterized by the shear modulus $\kappa_{\text{S}}$ while resistance to area dilatation is characterized by $C\kappa_{\text{S}}$ with $C$ much larger than unity for nearly area-incompressible membranes. For bending resistance the Helfrich model is used [90, 78]

[TABLE]

with $\kappa_{\text{B}}$ being the bending modulus, $H$ denoting the local mean curvature, $H_{0}$ the spontaneous curvature, $\kappa_{K}$ the Gaussian modulus and $K$ the Gaussian curvature. On the one hand, from a given energy functional, in-plane surface stresses and moments can be deduced by functional derivatives with respect to the strain tensor [74, 87] and curvature tensor, respectively [91, 92, 78]. On the other hand, in numerical algorithms the explicit introduction of stresses is typically bypassed and forces on the nodes discretizing the membrane are often computed directly by deriving a discretized version of the energy functional, here equations (7) and (8), with respect to the node positions [93, 50, 78]. This approach is used here as well for the passive elastic forces. For bending force calculation we use the method denoted by B in ref. [52].

For actively generated forces, however, this approach is not applicable since, due to their inherent non-equilibrium character, an energy functional cannot be defined in a strict sense. Instead, active contributions are usually given in terms of active stresses whose strength and direction can either be temporally constant or depend on additional quantities such a local concentration of motor proteins [6]. We will in the following construct a numerical method which explicitly applies these active stresses and derives the corresponding active forces on the discretized membrane nodes. These forces are then used to introduce a two-way coupling between active membrane dynamics and a surrounding hydrodynamic flow. For simplicity of demonstration and in order to connect to the analytical axisymmetric solutions of [3], we restrict ourselves to temporally constant active in-plane surface stresses. Time-dependent active stresses or active moments can be included without substantial modification of our numerical framework.

3.2 Force balance for deformable active membranes embedded in a 3D fluid

As sketched in figure 2 we consider a membrane immersed in an external fluid and enclosing an internal fluid, the cytosol. Coupling of the membrane to the internal and external fluid with the stress tensor $\sigma_{ij}^{in}$ and $\sigma_{ij}^{out}$, respectively, is described by the traction jump [76, 51, 88]

[TABLE]

with Latin subscripts denoting 3D coordinates and the components $n_{i}$ of the local unit normal vector onto the membrane pointing towards the external fluid. Transforming to the in-plane coordinate system the traction jump can - as a vector on the membrane in general - be decomposed into tangential and normal component

[TABLE]

Neglecting membrane inertia the traction jump between internal and external fluid is balanced by membrane forces (per unit area). In the present situation, these arise from elastic and active contributions. Looking ahead, the numerical method which we will construct in the following section will compute elastic forces $f_{\text{e}}^{\alpha}$ and $f^{n}_{\text{e}}$ via the classical discretized energy functional route, bypassing for simplicity the introduction of elastic stresses, while active forces need to be computed explicitly from active stresses and moments. It is thus convenient to write the force balance for a thin shell [73, 3, 87, 88, 6] as

[TABLE]

with active contributions to the in-plane and normal surface stress $t_{\text{a}}^{\alpha\beta},t_{n\text{a}}^{\alpha}$ and active contributions to the moments $m_{\text{a}}^{\alpha\beta},m_{n\text{a}}^{\alpha}$. The geometrical quantities are the curvature tensor $C_{\alpha\beta}^{\prime}$, the Levi-Civita tensor $\epsilon^{\prime}_{\alpha\beta}$ and the co-variant derivative $\nabla^{\prime}_{\alpha}$ taken on the deformed membrane as defined in the appendix. On the basis of equation (11) and (12) it becomes clear that the negative traction jump is the force exerted from the membrane on the fluid. We now consider the active contributions to the traction jump in the case of vanishing active moments ($m_{\text{a}}^{\alpha\beta}=m_{n\text{a}}^{\alpha}=0$) and vanishing active transverse shear surface stress ($t_{n\text{a}}^{\alpha}=0$)

[TABLE]

where we have used the definition of the co-variant derivative on the membrane (see equation (A1.20)) in the first line. To simulate the temporal dynamics and coupling to the external fluid of our discretized active membrane, we need to compute the forces on each membrane node corresponding to active in-plane surface stresses. According to equations (15) and (16) the curvature tensor, the Christoffel symbols, and the active in-plane surface stresses together with their derivatives have to be known locally on each node on a deformed surface. In the next section we will develop an algorithm to compute these quantities numerically for the discretized thin shell embedded in a 3D environment.

4 Algorithm for active force calculation on arbitrarily shaped discretized membranes

Our algorithm starts from active in-plane surface stresses and computes the corresponding active forces on the discretized membrane. The key ingredient of the algorithm is the discrete computation of geometrical properties on the discretized, deformed membrane, such as the curvature tensor or Christoffel symbols. This is achieved by a parabolic fitting procedure and using the force balance equations (15) and (16) as described in subsection 4.2 based on a local coordinate system which we introduce in subsection 4.1. With these forces the active membrane dynamics can be bidirectionally coupled to a surrounding fluid flow which is computed separately using either an overdamped dynamics or a Lattice-Boltzmann method.

4.1 Local coordinate system

As sketched in figure 3 each membrane node $\bm{r}_{c}$ has a neighborhood consisting of $N$ nodes, which are labeled in an ordered fashion around the central node. The choice of the starting node is arbitrary initially, but has to be retained through the simulation. In figure A1 we provide evidence that the choice of the starting node does not affect simulation results. At the site of the central node a local coordinate system ($\bm{e}_{\xi},\bm{e}_{\eta},\bm{n}$) can be defined, where we denote the in-plane coordinate vectors ($\bm{e}_{\alpha}$ ($\alpha=\xi,\eta$) in the appendix) by $\bm{e}_{\xi}$ and $\bm{e}_{\eta}$. We determine the local unit normal vector $\bm{n}$ by the mean weight by angle of the normal vectors on the surrounding triangles [94]. The first in-plane vector $\bm{e}_{\xi}$ is calculated using the vector from the central node to the first neighbor $\bm{x}_{1}=\bm{r}_{1}-\bm{r}_{c}$. The first step of the Gram-Schmidt orthogonalization is applied and the vector is normalized

[TABLE]

The second in-plane coordinate vector is calculated by the cross-product

[TABLE]

Using this method we can assign to every node at every time step a unique coordinate system which is adapted to the local deformed membrane geometry. We denote coordinates along $\bm{e}_{\xi}$ and $\bm{e}_{\eta}$ by $\xi$ and $\eta$, respectively. In a figurative sense, we co-move with a cytoskeletal filament initially positioned at $\bm{r}_{c}$ pointing towards $\bm{x}_{1}$.

4.2 Parabolic fitting

To obtain local geometrical quantities such as the curvature tensor on the deformed membrane we perform a parabolic fitting procedure based on the local coordinate system derived in the previous section. An arbitrary point $\bar{\bm{r}}$ in the neighborhood of the central node $\bm{r}_{c}$ can be expressed by a Taylor expansion around $\bm{r}_{c}$ up to the second order. For the $i$-th component of the vector $\bar{\bm{r}}$ with $i=x,y,z$ we obtain

[TABLE]

with $\xi$, $\eta$ being the coordinates along $\bm{e}_{\xi}$ and $\bm{e}_{\eta}$, respectively. Using this expression we can apply a parabolic fitting procedure (see also [52, 51] where a similar procedure has been used to compute passive bending forces) considering all $N$ neighboring nodes with a squared deviation from the fitted surface

[TABLE]

with $i=x,y,z$. By minimizing the $\chi_{i}^{2}$ we obtain the coefficients $A_{i}-E_{i}$. Using equation (19) we are able to analytically calculate the derivative of $\chi_{i}^{2}$ with respect to the coefficients $A_{i}-E_{i}$. The derivative of $\chi_{i}^{2}$ being zero in case of minimization then leads to a linear equation system for the $A_{i}-E_{i}$. This linear equation system is solved numerically in the simulation using LU decomposition. The paraboloid fitted to the neighborhood around $\bm{r}_{c}$ provides a good approximation of the local curvature for typical cell shapes [52]. By construction, the fitting coefficients equal the derivatives of the membrane parametrization vector $\bar{\bm{r}}$ with respect to local coordinates at the site of the central node

[TABLE]

or, equivalently,

[TABLE]

Thus, we are able to calculate geometrical quantities, as defined in the appendix, at the site of the central node in local coordinates with ${\tilde{\alpha}},{\tilde{\beta}},{\tilde{\gamma}}=\xi,\eta$, such as the metric tensor

[TABLE]

the curvature tensor

[TABLE]

and the derivatives of the metric tensor, e.g.,

[TABLE]

which are necessary for the calculation of the Christoffel symbols. The tensor $g^{{\tilde{\alpha}}{\tilde{\beta}}}$ is obtained by inverting the metric $g_{{\tilde{\alpha}}{\tilde{\beta}}}$ since $g_{{\tilde{\alpha}}{\tilde{\gamma}}}g^{{\tilde{\gamma}}{\tilde{\beta}}}=\delta_{{\tilde{\alpha}}}^{{\tilde{\beta}}}$. Thus, we can further calculate $C_{{\tilde{\alpha}}}^{\leavevmode\nobreak\ {\tilde{\beta}}}=C_{{\tilde{\alpha}}{\tilde{\gamma}}}g^{{\tilde{\gamma}}{\tilde{\beta}}}$ and $C^{{\tilde{\alpha}}{\tilde{\beta}}}=g^{{\tilde{\alpha}}{\tilde{\gamma}}}C_{{\tilde{\gamma}}}^{\leavevmode\nobreak\ {\tilde{\beta}}}$. We note that this procedure can be used without any restriction on the deformed membrane as well. Correspondingly, we obtain the metric tensor $g^{\prime}_{{\tilde{\alpha}}{\tilde{\beta}}}$, the curvature tensor $C^{\prime}_{{\tilde{\alpha}}{\tilde{\beta}}}$, and the Christoffel symbols on the deformed surface.

Equation (15) for the in-plane, active force includes a partial derivative $t_{\text{a}\leavevmode\nobreak\ ,\alpha}^{\alpha\beta}$, which is calculated using another parabolic fitting procedure. As done in equation (19) for the position in local coordinates we can expand the active in-plane surface stress $t^{{\alpha\beta}}$ around the components of the central node $t^{c{\alpha\beta}}_{\text{a}}$

[TABLE]

Corresponding to equation (20) we consider here the squared deviation from the expanded active stress $\bar{t}_{\text{a}}^{\tilde{\alpha}\tilde{\beta}}$

[TABLE]

with $t_{\text{a}}^{\tilde{\alpha}\tilde{\beta}\nu}$ being one component of the active in-plane surface stress of the $\nu$-th neighboring node in local coordinates. Thus, the corresponding fitting coefficients $A_{\text{a}}$ and $B_{\text{a}}$ represent the partial derivative of the active in-plane surface stress component with respect to the coordinate $\xi$ and $\eta$, respectively, whereas the coefficients $C_{\text{a}}\--E_{\text{a}}$ represent the second derivatives. We note that the matrix for the fit is the same as for the position. We perform one fitting procedure for each component, i.e., $\tilde{\alpha}=\xi,\eta$ and $\tilde{\beta}=\xi,\eta$.

With the method proposed here, we are able to calculate the co-variant derivative, which consists of a partial derivative and Christoffel symbols, and perform the contraction of active in-plane surface stress tensor and curvature tensor in the local coordinate system to determine the traction jump given by equations (15) and (16), which we repeat here for the deformed membrane in the local coordinate system

[TABLE]

in case of vanishing active moments and vanishing transverse shear stress. We note that the in-plane surface stress is the analogue of the Cauchy stress tensor in 3D [3, 88] and thus acting on the deformed membrane [95]. Equation (28) and (29) give the discrete forces (per area) on one membrane node. These forces then enter the fluid solver as illustrated in section 2.2 or are used for the relaxation dynamics in the overdamped limit (see appendix A3).

To convert the traction jump from equations (28) and (29) into a nodal force we use Meyer’s mixed area approach [78, 94]: in the case of non-obtuse triangles the area is calculated using Voronoi area $A_{\text{Voronoi}}$ defined by

[TABLE]

with the angles $\alpha$ and $\beta$ opposite to the edge connecting the central node and the $\nu$-th neighbor node within the adjacent triangles and in the case of obtuse triangles the midpoint of $\bm{x}_{\nu}$ that is opposite of the obtuse angle is chosen instead of the circumcenter point for each triangle.

4.3 Specification of active in-plane surface stress

Active stresses in membranes are often generated by ATP-fueled molecular motors "walking" along cross-linked rod-like structures such as actin filaments or microtubules. As a direct consequence, active stresses are often anisotropic with the direction of contractile/extensile stresses specified in a material frame moving and deforming along with the membrane itself. This naturally leads to a convenient specification of the active in-plane surface stress tensor $t_{\text{a}}^{\alpha\beta}$ in the local coordinate system introduced in section 4.1. Since the labeling of the neighbors around each node remains unchanged throughout the simulation, the distance vector $\bm{x}_{1}=\bm{r}_{1}-\bm{r}_{c}$ represents a material vector. Its normalized in-plane counterpart $\bm{e}_{\xi}$, given in equation (17), together with $\bm{e}_{\eta}$, given in equation (18), constitute the associated material frame. Imagining one cytoskeletal filament for an illustrative picture, the anchoring position of the filament is tracked by the node position $\bm{r}_{c}$, while the orientation of the filament is tracked by the fixed choice of the first neighbor $\bm{r}_{1}$. We again note that the choice of the starting node for the labeling, here 1, does not affect simulation results, as shown in figure A1.

The active in-plane surface stress tensor $t_{\text{a}}^{\alpha\beta}$ itself is not computed by our method but needs to be specified as an input quantity according to a corresponding constitutive law for the surface stress [3, 6]. Our algorithm allows for an arbitrary choice of active in-plane surface stress, subject to condition (14) for vanishing active moments, including spatially heterogeneous, anisotropic or time-dependent stresses. For the latter, coupling to a convection-diffusion model of active substances such as ATP or myosin (with the magnitude of $t_{\text{a}}^{\alpha\beta}$ proportional to local ATP/myosin concentration) is methodologically possible. Thus, if the concentration field of ATP/myosin is solved/prescribed on the membrane, e.g. by discretizing the convection-diffusion equation, the active stress can be calculated from the local concentration. The corresponding active forces and their coupling to the surrounding fluid dynamics are then straightforwardly achieved by the present algorithm, which allows for a spatially varying active stress. In the present work, we consider only temporally constant active stresses. Active stresses can thus conveniently be specified in the initial configuration of the membrane. For this, we first choose an intuitive coordinate system ($\bm{e}_{1}$,$\bm{e}_{2}$) appropriate for the initial shape of the undeformed cell membrane. The active in-plane surface stresses in this coordinate system are of the form

[TABLE]

where the mixed form with upper and lower index is chosen such that the right-hand side contains physical material-specific constants with dimensions of N/m [3]. Because of the dynamical deformation of the membrane, the active in-plane surface stress of equation (31) needs to be mapped into the local coordinate system at each node. This can be achieved by mapping the coordinate system ($\bm{e}_{1}$,$\bm{e}_{2}$) to the local coordinate system ($\bm{e}_{\xi}$,$\bm{e}_{\eta}$).

In many situations, the intuitive coordinate system ($\bm{e}_{1}$,$\bm{e}_{2}$) will correspond to cylindrical coordinates or spherical coordinates, e.g., for a rounded cell during mitosis [14]. By construction of the intuitive in-plane coordinate system and the local coordinate system both normal vectors, e.g., $\bm{e}_{r}$ for a cylinder or a sphere and the local $\bm{n}$ are equal initially. Thus, the initial in-plane coordinate system ($\bm{e}_{1}$,$\bm{e}_{2}$) can be converted directly into the local coordinate system ($\bm{e}_{\xi}$,$\bm{e}_{\eta}$) on every node individually

[TABLE]

with $\mathfrak{D}$ being a rotation matrix around the local unit normal vector $\bm{n}$, with $\alpha=1,2$ and ${\tilde{\alpha}}=\xi,\eta$. We thus obtain the active in-plane surface stress along the local coordinate vectors $\bm{e}_{\xi}$ and $\bm{e}_{\eta}$. The rotation in equation (32) is performed once at the beginning of a simulation. Since the local coordinate system $\bm{e}_{\xi}$, $\bm{e}_{\eta}$ is co-moving with the membrane material, the active in-plane surface stress tensor $t_{\text{a}}^{{\tilde{\alpha}}{\tilde{\beta}}}$, expressed in these coordinates, does not change over time.

The parabolic fitting procedure in equation (27) requires the difference in active stress between neighboring nodes. The differences in active stress are calculated in the intuitive coordinate system, in which the active stresses are prescribed. Because active stresses do not change in time, calculation of the differences can be done also once at the beginning of the simulation. Projection into the local coordinate system for each node is done in the same way as for active stresses in equation (32).

Along a certain direction, the cytoskeletal filaments may tend to contract or to expand, respectively, resulting in a contractile or extensile active in-plane surface stress [37]. The contractile or extensile nature of the cytoskeletal filaments manifests itself in the sign of the active stress, namely a positive (negative) stress corresponds to a contractile (extensile) nature. In active matter consisting of cytoskeletal filaments, active stresses often possess different signs in different directions. Imagining two anti-parallel polar filaments that are cross-linked by a motor protein walking in opposite directions on both filaments, a relative extensile shift of both filaments together with a lateral contraction occurs.

5 Validation

In this section we provide an in-depth validation of our algorithm to compute the dynamics of active membranes embedded in a 3D fluid. For this, we first compute the deformation obtained when an initially unstressed cylindrical membrane is subjected to a localized perturbation due to active stresses. Our numerical results are in excellent agreement with analytical predictions by Berthoumieux et al. [3] which were obtained using a Green’s function approach in the limit of small deformations. Next, we apply homogeneous active stresses again to an initially cylindrical membrane. In agreement with the analytical predictions of [3], we observe two kinds of axisymmetric instabilities. Going beyond the axisymmetric calculations of [3], our numerical method then predicts the existence of a third non-axisymmetric instability which, to the best of our knowledge, has not been observed thus far.

To account for the dynamics of the surrounding fluid, we employ two qualitatively different approaches. In the first approach, we use simple overdamped dynamics such that the surrounding fluid acts purely as a viscous frictional damping, see appendix A3. In the second approach, we consider the full fluid dynamics of the surrounding liquid by solving the Navier-Stokes equations using a Lattice-Boltzmann method, see section 2.1. Two-way coupling to the active membrane dynamics is provided by the Immersed-Boundary Method, as detailed in section 2.2.

5.1 Green’s function formalism

We here briefly recall the central analytical results of [3] which will be used to validate our numerical computations in the subsequent paragraphs. In [3] the active in-plane surface stresses have the following form

[TABLE]

on the initially unperturbed cylinder surface where the local coordinates $\alpha$ and $\beta$ correspond to polar coordinates $z$ and $\phi$, respectively (cf. 4.3). We note that the component $t_{\text{a}z}^{\leavevmode\nobreak\ z}$ being positive represents a contractile stress along the cylinder axis and $t_{\text{a}z}^{\leavevmode\nobreak\ z}$ being negative an extensile stress, as seen on the basis of the buckling instability reported for negative stress in ref. [3]. A positive $t_{\text{a}\phi}^{\leavevmode\nobreak\ \phi}$ represents a contractile stress in azimuthal direction, which causes a contraction of the cylinder. The latter becomes clear by considering the deformation caused by $T_{\text{a}}^{\phi}$ according to the Green’s function. In equation (33) $T_{\text{a}}$ represents an isotropic homogeneous background active stress while $T_{\text{a}}^{z}$ and $T_{\text{a}}^{\phi}$ are the amplitudes of point active stresses along each of the two coordinate axes. $\delta(z)$ is the Dirac delta distribution.

As in ref. [3], we focus on the radial deformation $u_{r}(z)$ of a cylinder with initial radius $R$ resulting from an azimuthal in-plane surface stress $T_{\text{a}}^{\phi}$. At the end of section 5.3 we perform a validation for an axial in-plane surface stress $T_{\text{a}}^{z}$. Berthoumieux et al. [3] consider a 3D elastic material and perform a projection onto the membrane resulting in the surface stretching modulus $S$ and the bending modulus $B$ as surface elastic parameters, together with the 3D Poisson ratio $\nu$. The strength of the homogeneous active in-plane stress is measured by the dimensionless number $g=\frac{T_{\text{a}}}{S}$ and bending forces are quantified by the relative bending modulus $b=\frac{B}{SR^{2}}$. The radial deformation of a cylindrical shell with radius $R$ is then given by the Green’s function $G_{r\phi}(z)$ as

[TABLE]

An expression for $G(z)$ is given in Fourier space [3] by

[TABLE]

with $G(z)$ being the inverse Fourier transform of $RG(q)$.

In our numerical method, elastic forces are derived from the Skalak and Helfrich energy functionals as defined in equation (7) and (8), respectively, which represent an accurate and widely used description of the elastic properties of biological cell membranes. In the limit of small deformations we show in section A2 that the Skalak and Helfrich model, which we use in this study, and the elastic model used by Berthoumieux et al. [3] are related by $S=\frac{2}{3}\kappa_{\text{S}}$ and $B=\frac{1}{2}\kappa_{\text{B}}$. In the following these relations are used to calculate the Green’s function in equation (35) and the corresponding deformation in equation (34) for comparison to our simulations. We furthermore set the bulk Poisson ratio $\nu=\frac{1}{2}$, corresponding to a 3D incompressible material, and the Skalak parameter $C=1$. Then the Skalak law becomes equivalent to the Neo-Hookean membrane law which is built to model a membrane made of a 3D incompressible material [88]. For the shear modulus $S$ we show in appendix A2 that the shear related in-plane surface stresses determining the Green’s function agree between the Skalak law and the model used by ref. [3] in the limit of small deformation. Although we can relate the parameter of Helfrich model to the model used in Berthoumieux et al. [3] by $B=\frac{1}{2}\kappa_{\text{B}}$ as illustrated in the appendix, the Helfrich model alters the in-plane surface stresses, which in turn alter the Green’s function. However, since we consider the limit $B\rightarrow 0$ this has no effect for the Green’s function used for validation of our simulations.

5.2 Singular active perturbation

We start by simulations of overdamped dynamics of an initially stress-free cylindrical shell subjected to a local increase in active in-plane surface stress. The active in-plane surface stress from equation (33) simplifies to

[TABLE]

as sketched in figure 4 a). Simulations are carried out with $R=1$, $\kappa_{\text{S}}=1$, $C=1$, $\kappa_{\text{B}}=0.001$, and $T_{\text{a}}^{\phi}=-0.01$ in simulation units. We employ two different triangulations of the cylindrical shape which are shown in 4 b) and c). In the coarse mesh, the axial distance between rings of nodes is $\Delta z=0.2$ and the azimuthal spacing is about $\Delta\phi=0.2$ radians. The finer mesh uses an axial spacing of $\Delta z=0.1$ and an azimuthal spacing of about $\Delta\phi=0.15$ radians.

In figure 5 we compare the final shape of the shell as observed in 3D simulations for the two different resolutions and compare it to the prediction of the Green’s function in equation (34). The analytical Green’s function shows a peak in deformation of finite width centered at the site of active in-plane surface stress perturbation ($z=0$) that decays with increasing distance, reaches a shallow minimum at around $z\approx\pm 0.7$ and then approaches zero for $z\to\pm\infty$. Although for both resolutions the resulting amplitude of the peak in deformation from the simulation is close to the Green’s function and both show a similar shape, we observe a significant deviation of the 3D simulation results from the theoretical expectation. Especially, the simulations cannot reproduce correctly the predicted shallow minima next to the main peak. We note that the deviations for singular perturbation also appear using Lattice Boltzmann/Immersed Boundary method instead of overdamped dynamics (results not shown).

This observed deviation can be explained by the idealized singular nature of the active in-plane surface stress which is impossible to accurately reproduce on a discretized membrane shape. Along the cylinder axis, i.e., along $z$-direction, only one ring of nodes with $z=0$ is attributed with finite active in-plane surface stress. Correspondingly, only these nodes are subjected to active forces and cause the neighboring nodes to move due to the elastic nature of the cylinder. We note that the central nodes - with active in-plane surface stress - experience a significantly larger deformation than their direct, adjacent neighbor nodes, as seen in figure 5. This steep gradient in deformation resulting in locally very large curvature cannot be completely resolved by the parabolic fit. Thus, the procedure fails to resolve the deformation caused by a point perturbation, although obtained deformations are similar in shape and nearly match the amplitude of the Green’s function. Given that in real applications, all perturbations can be expected to be non-singular, we proceed to investigate the behavior of our algorithm for spatially smooth active stresses.

5.3 Localized smooth perturbation

As a prototypical smooth distribution, we choose a Gaussian distributed active in-plane surface stress of the form

[TABLE]

where $T_{\text{a}}^{\phi}$ is again a constant amplitude. For this distribution, the predicted deformation can be obtained by superposing a distribution of Green’s functions (equation (34)) leading to the deformation

[TABLE]

We use the Green’s function in Fourier space given in equation (35) and the Fourier transform of the Gaussian

[TABLE]

Using the convolution theorem in Fourier space and transforming back to real space leads to

[TABLE]

The integral can be solved analytically in the limit of small bending rigidity $b\ll 1$ (which corresponds well with the chosen simulation parameters) to obtain

[TABLE]

Alternatively, we can solve the integral in equation (38) numerically, which does not lead to any differences for small bending modulus (results not shown).

In figure 6 a) we compare the 3D simulation results using overdamped dynamics for parameters $R=1$, $\kappa_{\text{S}}=1$, $C=1$ , $\kappa_{\text{B}}=10^{-5}$, and $T_{\text{a}}^{\phi}=-0.01$ to the analytical prediction of equation (39). The Gaussian distribution of active in-plane surface stress leads to a much smoother and broader peak of deformation than the singular perturbation of the previous subsection. Our simulation results are now in very good agreement with the theoretical prediction.

To go one step further, with applications in mind such as an active membrane in a flowing liquid, we now replace the simple overdamped fluid dynamics with a full Navier-Stokes dynamics solved by the Lattice-Boltzmann method and coupled to the active membrane via the Immersed-Boundary method as described in section 2. In figure 6 b) we compare simulation results obtained by LBM/IBM to the theoretical Green’s function. Again, our simulations are in very good agreement with the analytical theory. In figure 6 b) we include three sets of active in-plane surface stress and shear modulus with constant $\kappa_{\text{B}}=0.00018$ which are chosen such that the ratio of active in-plane surface stress and shear modulus $g$ remains constant. Thus, the Green’s function predicts identical deformation in all three cases which is indeed observed in our simulations.

We now proceed to the validation for a perturbation in $z$-stress, i.e., the active in-plane surface stress takes the form

[TABLE]

with the constant amplitude $T_{\text{a}}^{z}$. Corresponding to equation (38) the deformation can be obtained by [3]

[TABLE]

with $G_{rz}(s)=\frac{\nu}{RS}G(z)$ and takes in analogy to equation (39) the form

[TABLE]

In figure 6  c) we compare LBM/IBM simulation results for three different perturbations and shear moduli for $R=1$, $C=1$, and $\kappa_{\text{B}}=0.00018$ to the theory. Our simulation results are in very good agreement with the theory and by comparing figure 6 b) and c) we observe half the maximum deformation for $|T_{\text{a}}^{\phi}|=|T_{\text{a}}^{z}|$ which is indeed predicted by the theory for $\nu=\frac{1}{2}$.

As a further test of algorithm accuracy we now perform a convergence study based on the setup in figure 6 b) as well as in c). We compare the deformation obtained by simulation $u_{r}^{\text{sim}}$ and by Green’s function $u_{r}^{\text{Green}}$ by calculating the relative error per node defined as

[TABLE]

where $N_{z}$ denotes the number of nodes along the cylinder axis and the difference is in relation to the maximal deformation given by the Green’s function. In figure 6 d) we show the relative error per node in dependency of the number of membrane nodes. With increasing resolution the error per node steadily decreases in both cases.

The converge rate, however, is different for $T_{\text{a}}^{\phi}$ which decreases with slope $N_{z}^{-1}$ and $T_{\text{a}}^{z}$ which decreases more quickly with $N_{z}^{-2}$.

This may be due to the fact that the Green’s function is derived from linearized equations of motion [3], whereas simulations are also valid for larger deformations. A perturbation in $T_{\text{a}}^{z}$ shows half the maximal deformation as a perturbation in $T_{\text{a}}^{\phi}$ as predicted in equation (42) and thus is in better agreement with the Green’s function. The steady decrease in error demonstrates the accuracy of the presented algorithm for active force calculation.

From figure 6 we conclude that our algorithm presented in sec. 4 together with elastic force calculations gives reliable results for a reasonably smooth distribution of active in-plane surface stress. The very good accuracy of the predictions is achieved for the simple overdamped dynamics as well as for the substantially more complex and flexible combination of an active membrane with the lattice-Boltzmann/Immersed-Boundary method.

5.4 Homogeneous perturbation: instability diagram

To provide another test of our algorithm we perform simulations of a cylindrical membrane which now is subjected to a homogeneous active in-plane surface stress

[TABLE]

with $T_{\text{a}}=\text{const.}$ This situation corresponds to a membrane with constant motor protein activity and isotropic cortex architecture. We note again that positive $T_{\text{a}}$ represents contractile and negative $T_{\text{a}}$ extensile stress in $z$- or $\phi$-direction. Although they did not explicitly compute the deformation for this situation, ref. [3] predicts two unstable regions in $g$-$b$-parameter space for which the Green’s function could be shown to diverge. The predicted instability thresholds serve us as a further validation of our simulation method.

By varying both the relative active in-plane surface stress $g=T_{\text{a}}/S=3T_{\text{a}}/(2\kappa_{S})$ and the relative bending modulus $b=B/(SR^{2})=3\kappa_{B}/(4\kappa_{S}R^{2})$ we obtain the phase diagram in figure 7 and compare it to the predicted instability thresholds given by Berthoumieux et al. [3]. On the right, for large $g$, Berthoumieux et al. [3] predict an instability occurring for $T_{\text{a}}>\frac{4}{3}\kappa_{S}(1-\nu^{2})$ shown by the vertical orange line in figure 7. Indeed, for simulations in this range we observe an instability with local contraction of the cylinder (see inset of figure 7 for a shape illustration). The threshold obtained by our simulations closely matches the analytically predicted threshold. This instability is analogous to a Rayleigh-Plateau instability of a liquid jet with the positive, contractile active in-plane surface stress playing the role of the surface tension. To the left of the threshold a fairly large region is observed in which the initial cylindrical shape remains stable. For negative $g$ (extensile stress), [3] predicts a buckling instability when $T_{\text{a}}<-2\sqrt{\frac{\kappa_{B}\kappa_{S}}{R^{2}}}$. To compare to our simulations, we prescribe an active stress only along the cylinder axis, i.e., $t_{\text{a}z}^{\leavevmode\nobreak\ z}=T_{\text{a}}$ and $t_{\text{a}\phi}^{\leavevmode\nobreak\ \phi}=0$. This corresponds to a contracting, cylindrical membrane and the resulting shape beyond the threshold is illustrated in figure 7 at the bottom left. Our simulations agree very well with the predicted instability onset depending on the relative active in-plane surface stress and the relative bending modulus.

In addition, we investigate the transition to buckling in more detail and carry out simulations imposing an active tension in azimuthal and axial direction $t_{\text{a}z}^{\leavevmode\nobreak\ z}=t_{\text{a}\phi}^{\leavevmode\nobreak\ \phi}=T_{a}$, which corresponds exactly to the scenario considered by [3]. In figure 8 we compare these simulations (top row) with the ones in figure 7 (bottom row), respectively. At large (negative) active stress in a) we observe an instability in both simulations, however only the simulation with a purely axial stress clearly corresponds to a buckling instability. The instability in the top row exhibits a non-axisymmetric character. For slightly smaller active stress in 8 b) the non-axisymmetric instability remains for the isotropic stress, but the membrane becomes stable for $z$-stress only. Decreasing the active stress further, we observe a stable cylindrical membrane in both cases, as shown in c).

From figure 8 we conclude that an additional instability is present (not related to buckling) caused by a finite azimuthal stress $t_{\text{a}\phi}^{\leavevmode\nobreak\ \phi}$. This additional instability induces non-axisymmetric deformations and thus is not observed in the axisymmetric treatment of ref. [3] nor in our axisymmetric simulations for which also the onset of the buckling instability is in exact agreement with the analytical prediction as shown in figure S2 of the Supporting Information.

So far, we have validated our algorithm to agree with theoretical predictions on the basis of the Green’s function. To obtain a validation in the non-linear regime beyond the Green’s function, we compare our 3D simulations to simulations of an axisymmetric membrane, as detailed in the Supplemental Material [96]. We do this by comparing the membrane shape in the case of a buckling instability. In figure 9 we show the deformation obtained from the axisymmetric simulation and compare it to three simulations using 3D LBM/IBM. The three 3D simulations are done for different resolutions $\Delta z/R$. We use the non-dimensional parameters $g=-0.75$ and $b=0.01125$. On the one hand all three 3D simulations show the same deformation and wavelength and on the other hand they agree in the wavelength with the axisymmetric simulation method. The wavelength of the buckling instability is about 1.75$R$ for both very different simulation techniques. In case of small bending elasticity Berthoumieux et al. [3] predict a wavelength at the threshold, where the denominator of the Green’s function in eq. (35) becomes zero for finite wave vector $q$, of $\lambda=1.9R$ which is reasonably close to the value observed with simulations. The difference may arise from the periodicity of the shell in our simulations and/or from finite bending together with being beyond the threshold. Nevertheless, the excellent agreement between the axisymmetric and 3D simulations provides strong evidence for the reliability of our algorithm also in the range of large deformations.

6 Model application: cell division in suspending fluid

In the following we present a first application of our method including fluid flow. This illustrates the versatility and applicability of our combined LBM/IBM method for active cell membranes. For this, we consider a dividing ellipsoidal cell. Except the fact that we employ an elastic rather than a viscous cortex, our setup resembles the situation of cell cytokinesis [16, 17, 19, 20, 21]. Cell cytokinesis as part of cell division is a prominent subject of active matter research in biological physics [13, 18, 97]. Most previous studies, e.g. [20] or [21], investigated the dynamics of cell cytokinesis for an axisymmetric membrane without considering internal fluid flow. Refs. [98] and [99] analyzed the flow field inside a dividing cell where the contractile ring is modeled as an additional force using the immersed boundary method and phase field model, respectively. Ref. [100] analyzed the flow field by means of the phase field model as well, but considered the cortical ring as shrinking elastic loop. Here, we consider cell division triggered by active stresses including an external flow leading to a fully 3D asymmetric membrane shape during the division. We first analyze the flow field dynamically evolving inside the dividing cell surrounded by a quiescent medium in section 6.1 and then extend this setup by considering a dividing cell in an external shear flow in section 6.2.

6.1 Flow field inside a dividing cell

We consider a prolate ellipsoidal cell of diameter $7\leavevmode\nobreak\ \mu$m and length $14\leavevmode\nobreak\ \mu$m which is endowed with an isotropic active stress $T_{\text{a}}=8\times 10^{-5}\leavevmode\nobreak\ \text{N/m}$. In addition, in an interval of $\Delta\theta\approx\frac{\pi}{12}$ around the equator the active stress in azimuthal direction is increased by a factor of six according to a step function. The membrane is endowed with shear elasticity $\kappa_{S}=5\cdot 10^{6}$ N/m, $C=1$, and bending elasticity $\kappa_{B}=2\cdot 10^{-19}$ Nm which are in the range of typical cell membranes.

We present the initial 3D membrane shape in figure 10 a) and a deformed 3D shape in figure 10 b). In both cases we illustrate the region with increased active stress along $\phi$-direction by the red shaded area. The shape and the developing flow field in the central plane, which includes the long axis of the ellipsoid, are shown over time in figure 10 c) to h).

Here, the active stress triggers active deformation of the membrane which in turn triggers fluid flow inside the cell. We note that the fluid velocity at the position of the membrane corresponds to membrane motion, which moves with local fluid velocity due to no-slip condition as described in section 2.2.

At the beginning the membrane contracts around the poles (left and right in figure 10 c)) due to the isotropic contractile active stress. The contraction at the poles causes a rounding which triggers a flow field pointing away from the poles at the beginning. Simultaneously, the membrane starts contracting around the equator, see figure 10 c) and d). In figure 10 d) four vortices are present around the equator and four at the corners of the figure. After some time the contraction at the poles stops, see figure 10 e) and f) and only the contraction at the equator proceeds. With progressive contraction a flow away from the midplane towards the poles develops. As it is visible in figure 10 g) and h), the site of maximal velocity towards the poles is located at $x\approx\pm 10$. Thus, it does not coincide with the center of the two spheroids pinching off at $x\approx\pm 18$, but is rather shifted towards the equatorial plane at $x=0$. Approaching the poles of the ellipsoid the velocity decreases.

6.2 Dividing cell in shear flow

In the previous section we considered an initially quiescent fluid. Here, we go one step further and include an externally driven flow interacting with the membrane. We apply a shear flow with a shear rate $\dot{\gamma}\approx 1400$ s*-1*. All other parameters are the same as in the previous section.

The cell membrane as illustrated in figure 11 a) deforms now due to the local increase in active stress (red shaded area in figure a) and b)) but also due to the external shear flow. This becomes visible by the non-symmetrically deformed membrane in figure 11 b). Figures 11 c) to h) show the flow field relative to the shear flow, i.e., from each velocity vector the corresponding background flow is subtracted. The time evolution of the cell shape and the flow field is also illustrated in the Supplemental Video 1. In contrast to the previous section, the shear flow triggers an asymmetric deformation and in turn an asymmetric flow inside the cell. The rounding together with the flow from the poles towards the equator is less pronounced (compare figure 11 d) to figure 10 d)).

The flow field inside a dividing cell suspended in a shear flow shows how the actively deforming membrane couples to a background flow and imposes perturbation on the shear flow. Both the active stress present in the cell cortex as well as the external flow trigger membrane deformation.

7 Conclusion

We presented a computational algorithm to compute the dynamical deformation of arbitrarily shaped active biological (cell) membranes embedded in a 3D fluid. Active stresses in cells typically arise from the activity of motor proteins. Constitutive equations for active stresses in membranes have been developed recently [30, 6] in the framework of differential geometry and form the theoretical basis for our computational method. The membranes are discretized by a set of nodes connected via flat triangles. The key ingredient of our algorithm is the computation of the active force acting on each node starting from prescribed active stresses via a parabolic fitting procedure on the deformed membrane. Besides active forces, the method also includes passive elastic forces derived from the well-established Skalak and Helfrich models for cell membranes. In simple cases, the surrounding fluid can be considered a purely frictional medium, such that the membrane nodes follow simple overdamped dynamics in time. For more realistic situations, we introduced a powerful and versatile coupling between the active membrane and the surrounding fluid via the Immersed-Boundary method. This technique incorporates the full Navier-Stokes dynamics for the surrounding liquid, solved here via the Lattice-Boltzmann method. Thus, our method allows us to go beyond the determination of equilibrium shapes of active elastic membranes and allows for simulations of dynamically deforming biological cells immersed in an external flow.

We successfully validated our algorithm for two distinct situations of an elastic, initially cylindrical membrane: (i) a local, Gaussian distributed active stress and (ii) a homogeneous active stress. For (i) our numerical results are in excellent agreement with the analytically predicted deformation of Berthoumieux et al. [3] and show convergence with increasing resolution. Overdamped dynamics and IBM/LBM dynamics are in good agreement. For (ii) we recovered both the buckling as well as the Rayleigh-Plateau-like instability predicted by [3]. Comparison to our own numerical solutions of the axisymmetric problem shows very good agreement with the full 3D algorithm also in the regime of large deformation not covered by [3]. In addition, our computations reveal the existence of a thus far unobserved non-axisymmetric instability in the case of extensile axial and azimuthal stresses. In order to illustrate the versatility of our method, we analyzed the flow field inside an elastic, dividing cell membrane in shear flow. This represents the first investigation of the dynamic two-way coupling between active deformations and externally driven fluid flow. In this work, we considered temporally constant active stresses, but the inclusion of time-dependent active stresses computed, e.g., by a convection-diffusion model of active substances, is straightforwardly possible.

Our computational method significantly extends the range of physical problems to which existing active membranes theories can be applied. First, it is not restricted to simple shapes such as cylinders or spheres (or small deformations thereof) allowing efficient and accurate treatment of arbitrary membrane shapes and deformations. Second, our method couples the active membrane dynamics to the full Navier-Stokes dynamics of the surrounding fluid.

In the described LBM/IBM scheme a viscosity contrast of the inner and outer fluid - as it is well known in case of red blood cells [50] - can furthermore be incorporated, which allows for a even more realistic model of living cells.

This opens up a wide range of applications in external flows such as active cells in the blood stream or active cellular compartments in cytoplasmic streaming flows which currently remain largely unexplored. A particularly interesting application could be the formation of platelets from megakaryocytes which, according to a set of recent experiments [101, 102], crucially depends on the interplay between active processes and external flows.

8 Acknowledgments

C. B. thanks the Studienstiftung des deutschen Volkes for financial support and acknowledges support by the study program "Biological Physics" of the Elite Network of Bavaria. This project is part of the collaborative research centre TRR225 (subproject B07) funded by the German Research Foundation (DFG). We gratefully acknowledge computing time provided by the SuperMUC system of the Leibniz Rechenzentrum, Garching, as well as by the Bavarian Polymer Institute and financial support from the Volkswagen Foundation.

Appendix

A1 Membrane in thin shell formulation

In this appendix, we summarize the necessary basics and conventions of differential geometry on thin shells used in this work. For a more detailed description we refer the reader to refs. [73, 84]. The 2D manifold in general is parametrized by two coordinates $s^{1},s^{2}$. We denote vector components on the manifold by Greek letters $\alpha,\beta=1,2$ and vector components in Euclidean space by Latin letters $i,j,k=x,y,z$. Moreover, we use the Einstein summation convention, i.e., double occurrence of an index in sub and superscript implies a sum over this index. A partial derivative with respect to $s^{\beta}$ is denoted by a comma, i.e., for an arbitrary vector $\bm{v}$ this implies $v_{\alpha,\beta}=\partial v_{\alpha}/\partial s^{\beta}$.

The membrane in the undeformed state is parametrized by the vector

[TABLE]

From the local in-plane coordinates

[TABLE]

the metric tensor is defined by

[TABLE]

with $g=\det(g_{\alpha\beta})$. The inverse metric $g^{\alpha\beta}$ can be obtained by the expression

[TABLE]

with $\delta_{\alpha}^{\beta}$ denoting the Kronecker delta being 1 for $\alpha=\beta$ and 0 otherwise. For a general vector $v_{\alpha}$, $v^{\alpha}$ or tensor $t_{\alpha\beta}$, $t_{\alpha}^{\leavevmode\nobreak\ \beta}$, $t^{\alpha}_{\leavevmode\nobreak\ \beta}$, $t^{\alpha\beta}$ subscript indices denote co-variant components and superscript indices denote contra-variant components. An index can be raised by

[TABLE]

or lowered by

[TABLE]

A surface element is defined by $dS=\sqrt{g}ds^{1}ds^{2}$. The Christoffel symbols are given by

[TABLE]

The in-plane coordinate vectors $\bm{e}_{\alpha}$ provide a local coordinate system together with the unit normal vector on the membrane

[TABLE]

Considering the local unit normal vector allows for the definition of the curvature tensor

[TABLE]

A co-variant derivative of an arbitrary vector $v^{\alpha}$ or tensor $t^{\alpha\beta}$ defined on the membrane is given by

[TABLE]

On the membrane the Levi-Civita tensor is given by

[TABLE]

An arbitrary vector $\bm{v}^{\alpha}$ defined on the manifold with respect to $\bm{e}_{\alpha}$ with $\alpha=1,2$ can be decomposed into an in-plane and a normal contribution

[TABLE]

with $t^{\alpha\beta}$ being the component of a tensor and $t_{n}^{\alpha}$ the component of a vector.

A connection from the in-plane coordinates to the Euclidean coordinates can be drawn by the expression

[TABLE]

with $e_{\alpha}^{i}$ being the $i$-th component of $\bm{e}_{\alpha}$ and $\bm{E}_{i}$ being the $i$-th Euclidean unit vector. A three dimensional tensor $t_{ij}$ can be projected onto the membrane via

[TABLE]

External or internal forces may lead to a deformation of the membrane characterized by the deformation field $\bm{u}$. The membrane in the deformed state is parametrized by

[TABLE]

We denote all vectors or tensors that are evaluated on the deformed surface by a prime. Corresponding to the change of local coordinate vectors $\bm{e}_{\alpha}^{\prime}=\bm{X}_{,\alpha}^{\prime}$ and normal vector $\bm{n}^{\prime}$ both the metric tensor and the curvature tensor changes to

[TABLE]

The Christoffel symbols have to be computed using $g_{\alpha\beta}^{\prime}$ and the co-variant derivative becomes

[TABLE]

A2 Elastic in-plane surface stresses

In the following we compare the elastic in-plane surface stresses used in ref. [3] to those obtained for Skalak energy density in equation 7. We consider the displacement vector in equation A1.16 decomposed into axial and normal deformation

[TABLE]

For the given energy density in equation 7 the in-plane surface stresses are obtained by [74, 103, 104]

[TABLE]

with the invariants

[TABLE]

and with $J=\sqrt{I_{2}+1}$. Using equations A1.16, A1.17 and the deformation in equation A2.1 we obtain for the metric on the deformed membrane in the limit of small deformations

[TABLE]

and for the in-plane surface stresses

[TABLE]

These equations can be compared to the tensions in equation (11) of ref. [3] for the elastic model used by Berthoumieux et al. [3]. The latter is based on the Hooke’s law in 3D which is projected onto the membrane. By comparing the stresses we find agreement in the limit of small deformations for the relation of the stretching modulus $S$ to the shear modulus used in equation (7) of

[TABLE]

and $C=1$. This relation is used for the Green’s function to match both elastic models.

We furthermore can compare the bending modulus used in Berthoumieux et al. [3] to the one appearing in the Helfrich energy in equation (8) [78]. Berthoumieux et al. [3] defines $B=Eh^{3}/(24(1-\nu^{2}))$ with the Youngs modulus $E$ and the thin shell height $h$. Comparing this to the expression $\kappa_{B}=Eh^{3}/(12(1-\nu^{2}))$ of Pozrikidis [105] we obtain the relation

[TABLE]

A3 Ovedamped dynamics method and simulation analysis

A3.1 Overdamped dynamics

As a different approach than the LBM/IBM, we use a model program based on overdamped dynamics to solve for the final, equilibrium shape of the membrane in case of validation. The resulting active and elastic forces $\bm{F}$ calculated for every node enter the equation of motion of the corresponding node $\bm{r}_{c}$ which is given by

[TABLE]

where $\gamma$ is a friction coefficient. We solve the equations for all nodes using Euler integration scheme. We fix the nodes at the boundaries of the cylinder by harmonic springs of strength $1000\kappa_{S}$. This results in nearly vanishing deformation at the boundary of the cylinder.

A3.2 Simulation analysis

To obtain the shape shown for example in figure 6 a) or b) we average the final, radial deformation over all nodes at certain $z$-position. Due to the averaging and inherent errors in the bending algorithm [52] the deformation $u_{r}$ does not reach exactly zero far away from the perturbation in active stress, but show a constant offset of $5\times 10^{-4}$, which we eliminate in the figures.

A4 Choice of reference neighbor does not influence simulation results

In section 4.1 we explain how the local coordinate system is built using the first neighboring node as reference for the first in-plane coordinate vector $\bm{e}_{\xi}$. This fact is also important for the projection of the active stress into the local coordinate system, as mentioned in section 4.3. In the following we show that the choice of the reference neighboring node is arbitrary and does not influence simulation results.

Therefore, we first consider the cylindrical membrane subjected to a homogeneous perturbation in active in-plane surface stress along $z$-direction $T_{\text{a}}^{z}$. We use exactly the setup analyzed in figure 6 c) with $T_{\text{a}}^{z}=0.01$ and $\kappa_{S}=1.0$. We systematically change the reference neighbor node serving for construction of the local coordinate system as illustrated in figure A1 a). Figure A1 b) shows that the deformation obtained in simulations is the same for every choice of reference neighbor node and agrees very well with the theoretically obtained Green’s function. In case of $T_{\text{a}}^{z}$ perturbation the in-plane derivatives of the active in-plane surface stress are crucial and trigger the deformation. Thus, figure A1 provides evidence for the correct calculation of derivatives regardless of the choice of reference neighbor node and furthermore shows the choice is arbitrary.

Furthermore, in order to take into account a more complex membrane geometry, we use the setup in figure 10 and compare the dynamics for five different choices of reference neighboring node. Choosing the sixth neighbor is not possible for the ellipsoidal geometry as the surface tiling requires at least 12 nodes with five neighbors only. We consider an initially ellipsoidal cell membrane endowed with active stress, which in azimuthal direction increases around the equator. The increased active stress triggers the cell membrane to contract. In figure 10 we show the dynamically evolving flow field inside the cell. Here, we redo the simulation and for each simulation we choose a different neighbor node to build the local coordinate system, which moves with the deforming membrane in time. In figure A1 c) we show the five different choices.

In figure A1 d) we consider the point in time corresponding to figure 10 g) and show the radial position of all membrane nodes as function of the position along the axis. All simulations with different reference neighbor show the same membrane shape and are in very good agreement. We note that slight deviations occur at the site of strongest indentation, which results also in the strongest curvature. At this position, the parabolic fit is not capable of covering the strongly deformed membrane shape completely and thus slight deviations occur. In figure A1 e) we do the same for a dividing cell in shear flow corresponding to figure 11 g). We show the node positions within the plane containing the long axis of the cell and the shear axis of the external flow. All simulations with different reference neighbor show the same membrane shape in shear flow and are in very good agreement. Thus, figure A1 d) and e) provide evidence that the choice of reference neighbor is indeed arbitrary and does not alter the simulation in case of a dynamically deforming membrane coupled to a suspending fluid.