In general, the two dimensional thin shell representing the cell cortex is parametrized by the three-dimensional vector X (s ¹, s ²) which depends on the two surface coordinates s ¹, s ². The latter determine the in-plane position on the thin shell. From the parametrization X two in-plane coordinate vectors pointing locally along the thin shell are derived e 1 = ∂ X ∂ s 1 , e 2 = ∂ X ∂ s 2 , (1) and using those a unit normal vector on the thin shell that points outwards can be defined n = e 1 × e 2 | e 1 × e 2 | . (2)

The metric tensor of the thin shell is defined by g α β = e α ⋅ e β , (3) and its curvature tensor by C α β = − ∂ α ∂ β X ⋅ n , (4) with Greek indices α, β, γ = 1, 2 referring to the in-plane coordinates s ¹, s ². In the following, we use the Einstein convention where a double occurrence of an upper (contravariant) and lower (covariant) index implies a sum. A lower (covariant) index can be raised with the metric tensor, e.g., for a general vector component a ^β = a _α g ^αβ or a general tensor t _α ^β = t _αγ g ^γβ . Correspondingly, an upper (contravariant) index can be lowered by a _β = g _βα a ^α or t _αβ = t _α ^γ g _γβ . For a symmetric tensor t _αβ = t _βα the order of upper and lower indices becomes obsolete in mixed notation, i.e., it can be written as t α β .

We denote the partial derivative along the in-plane coordinates of a general function f by ∂ _α f and the covariant derivative by ∇ _α f. The covariant derivative is defined for a scalar by ∇ _α f = ∂ _α f, for a vector component by ∇ α a β = ∂ α a β + Γ α γ β a γ , and for a general tensor by ∇ α t β γ = ∂ α t β γ + Γ δ α β t δ γ + Γ α δ γ t β δ , where Γ α β γ denotes the Christoffel symbols of the second kind [54]. Later, we denote the contraction of the covariant derivative with the metric tensor in contravariant components as ∇ ^β = g ^βγ ∇ _γ . The covariant derivative of the in-plane and normal coordinate vector is given by the equations of Gauss and Weingarten ∇ α e β = − C α β n , ∇ α n = C α β e β , (5) respectively.

The key quantity of interest for the viscous active cortex is the cortical velocity field v . It is defined on the cortex, thus depending on (s ¹, s ²), but is a three-dimensional vector. Compared to Ref. [45] we do not use a Hodge decomposition of the velocity field, but treat it as a vector field expressed in Cartesian coordinates. Therefore, it can be decomposed in the local coordinate system on the thin shell by projection onto the in-plane and normal coordinate vector v α = v ⋅ e α a n d v n = v ⋅ n , (6) with its in-plane components v ^α and its normal component v ⁿ such that v = v α e α + v n n . (7)

The total in-plane vectorial component of the velocity field can further be determined by v t = 1 − n ⊗ n ⋅ v , (8) with the projector n ⊗ n onto the normal vector, which is the outer product of the normal vector with itself, and with 1 the unit matrix. The velocity gradient on the thin shell is defined as [40]. v α β = 1 2 ∇ α v β + ∇ β v α + C α β v n , (9) which is equivalent to the definition in vector notation v α β = 1 2 ∂ α v ⋅ e β + ∂ β v ⋅ e α = 1 2 ∇ α v ⋅ e β + ∇ β v ⋅ e α . (10)

Forces in the cortex, e.g., arising due to gradients in the velocity, are described by a stress tensor, in analogy to three-dimensional hydrodynamics [53, 56]. Here, the stress tensor is defined on the thin shell and therefore denoted as surface stress t ^α [42, 53] consisting of two vectors (α = 1, 2) and having dimensions 3 × 2. The surface stress in turn can be further decomposed into the in-plane surface stress t ^αβ , a 2 × 2-tensor, and the normal surface stress t n α [42] as t α = t α β e β + t n α n . (11)

For vanishing normal surface stress t n α = 0 , which is typically the case for an infinitely thin surface which does not sustain internal bending moments [40], the force balance for a thin shell in the presence of a pressure difference P between the cell’s interior, the cytoplasm, and its external environment, but in absence of other external forces, becomes. ∇ α t α β = 0 , (12) − C α β t α β = − P . (13)

The pressure difference P accounts for the incompressibility of the cytoplasm, which enters as an additional constraint as detailed in section 3.2.

We consider a viscous cortex with planar shear η _s and bulk viscosity η _b , which is subject to active forces described by the active surface stress tensor ζ ^αβ . Accordingly, the constitutive equation reads [40]. t α β = t v α β + ζ α β , (14) with the viscous surface stress tensor determined by t v α β = 2 η s v ̃ α β + η b v γ γ g α β , (15) where v ̃ α β = v α β − 1 2 v γ γ g α β . The expression in Eq. 15 is equivalent to the viscous stress tensor in three dimensions [56]. The active surface stress ζ ^αβ describes internal forces in the cortex which stem from active processes, such as motor proteins walking along cytoskeletal filaments by conversion of chemically stored energy [35, 37, 40]. The active surface stress ranges around 10⁻⁴N/m [57–61]. A typical value of the shear viscosity of the cell cortex is η _s = 27 × 10⁻⁴ Pa s m [62]. For the plasma membrane of a red blood cell a shear viscosity of about η _s ≈ 3 × 10⁻⁷ Pa s m has been reported [63, 64], which is about four orders of magnitude smaller than a typical cortex viscosity and therefore the cortical viscosity dominates. We expect in general the two viscosities to be of similar order of magnitude, if the cortex can be seen as a thin layer of homogeneous material.

In the following, we assume that the distribution of active surface stress ζ ^αβ in the cortex is known, but our method allows for a potential coupling of active surface stress to a concentration field, e.g., of myosin. We consider a coupling to a passive immersing fluid, where the pressure difference P balances the average active in-plane surface stress. Despite these contributions we consider the cortex to be free of other external forces or external torques. Inserting the constitutive law for in-plane surface stress tensor (14) into the force balance Eqs 12, 13 we obtain for the cortex in absence of any other external forces. 2 η s ∇ α v α β + η b − η s ∇ α v γ γ g α β = − ∇ α ζ α β , (16) − 2 η s C α β v α β − η b − η s v γ γ C α β g α β + P = C α β ζ α β , (17) where the Eq. 16 are the tangential force balance equations and Eq. 17 is the normal one. Our goal is to solve the force balance Eqs 16, 17 for the velocity field and pressure difference due to a given active surface stress distribution ζ ^αβ .

The aim of this section is a formulation of the force balance Eqs 16, 17 in terms of the velocity and its derivatives expressed in Cartesian vectors. This is the direct basis of our numerical algorithm, which will be detailed in section 3, and we obtain this by combining section 2.1 and section 2.2. Motivated by the parabolic fitting procedure to determine curvature and derivatives on a discrete thin shell [42, 65], we start with a parabolic expansion of the velocity vector in in-plane coordinates (s ¹, s ²) in the surrounding of a given position r _ν on the thin shell v s 1 , s 2 = v ν + ∇ 1 v ︸ A v s 1 + ∇ 2 v ︸ B v s 2 + 1 2 ∇ 1 ∇ 1 v ︸ C v s 1 s 1 + 1 2 ∇ 2 ∇ 2 v ︸ D v s 2 s 2 + ∇ 1 ∇ 2 v ︸ E v s 1 s 2 , (18) with v ν = v ( r ν ) = v s 1 = 0 , s 2 = 0 and A _v , B _v , …, E _v being the first and second derivative of the velocity field v with respect to in-plane coordinates at position r = r _ν . Using the parabolic expansion of the velocity field in Eq. 18, we can first directly express the velocity gradient (10) evaluated at r = r _ν in terms of the velocity derivatives A _v , B _v , …, E _v v 11 = A v ⋅ e 1 , v 22 = B v ⋅ e 2 , (19) v 12 = v 21 = 1 2 A v ⋅ e 2 + B v ⋅ e 1 . (20)

The trace of the velocity gradient is v γ γ = v α β g α β = g 11 v 11 + 2 g 12 v 12 + g 22 v 22 = g 11 A v ⋅ e 1 + g 12 A v ⋅ e 2 + B v ⋅ e 1 + g 22 B v ⋅ e 2 . (21)

Next, we can calculate the derivative of the velocity gradient tensor in Eq. 10 which, using the equations of Weingarten and Gauss (5), becomes ∇ α v β γ = 1 2 ∇ α ∇ β v ⋅ e γ − 1 2 C α γ ∇ β v ⋅ n + 1 2 ∇ α ∇ γ v ⋅ e β − 1 2 C α β ∇ γ v ⋅ n . (22)

Using that ∇ _α  g ^βγ = 0, we can calculate the derivative of the contravariant components of the velocity gradient by ∇ α v β γ = g β ϵ g δ γ ∇ α v ϵ δ . (23)

Furthermore, we calculate the derivative of the velocity gradient’s trace via ∇ α v γ γ = ∇ α g γ β v β γ = g γ β ∇ α v β γ = g 11 ∇ α v 11 + g 12 ∇ α v 12 + g 21 ∇ α v 21 + g 22 ∇ α v 22 using the formula for the gradient in Eq. 22. Writing down these equations for fixed indices α, β, γ = 1, 2, they can be explicitly written in terms of A _v , B _v , …, E _v .

In the final step, we aim for a formulation of the force balance equations Eq. 16, Eq. 17 in terms of the velocity derivatives in Cartesian coordinates. Rearranging equations Eq. 16, Eq. 17 using v β α = v β γ g γ α and g α β g β γ = δ α γ leads to 2 η s g γ α ∇ α v β γ + η b − η s ∇ α v γ γ δ β α = − ∇ α ζ β α (24) − 2 η s C α β v α β − η b − η s v γ γ C δ δ + P = C β α ζ α β . (25)

The components of the force balance equation in this form can be explicitly expanded using the expressions derived above. Writing down each component on its own and collecting terms with respect to A _v , B _v , …, E _v we end up with the first tangential force balance equation in the form − 2 η s g 11 C 11 n − 3 η s g 21 C 12 n − η s g 22 C 22 n − η b − η s g 11 C 11 n − η b − η s g 12 C 12 n ⋅ A v + − η s g 21 C 11 n − η s g 22 C 21 n − η b − η s g 12 C 11 n − η b − η s g 22 C 12 n ⋅ B v + + 2 η s g 11 e 1 + η s g 21 e 2 + η b − η s g 11 e 1 + η b − η s g 12 e 2 ⋅ C v + + η s g 22 e 1 ⋅ D v + + 3 η s g 21 e 1 + η s g 22 e 2 + η b − η s g 12 e 1 + η b − η s g 22 e 2 ⋅ E v = − ∇ 1 ζ 1 1 − ∇ 2 ζ 1 2 , (26) the second tangential force balance equation − η s g 11 C 12 n − η b − η s g 12 C 22 n − η s g 12 C 22 n − η b − η s g 11 C 21 n ⋅ A v + − η s g 11 C 11 n − 3 η s g 12 C 21 n − 2 η s g 22 C 22 n − η b − η s g 12 C 21 n − η b − η s g 22 C 22 n ⋅ B v + + η s g 11 e 2 ⋅ C v + + η s g 12 e 1 + 2 η s g 22 e 2 + η b − η s g 12 e 1 + η b − η s g 22 e 2 ⋅ D v + + η s g 11 e 1 + 3 η s g 12 e 2 + η b − η s g 11 e 1 + η b − η s g 12 e 2 ⋅ E v = − ∇ 1 ζ 2 1 − ∇ 2 ζ 2 2 , (27) and the normal force balance equation − 2 η s C 11 e 1 − 2 η s C 12 e 2 − C η b − η s g 11 e 1 − C η b − η s g 12 e 2 ⋅ A v + − 2 η s C 12 e 1 − 2 η s C 22 e 2 − C η b − η s g 12 e 1 − C η b − η s g 22 e 2 ⋅ B v + P = C 1 1 ζ 1 1 + C 1 2 ζ 2 1 + C 2 1 ζ 1 2 + C 2 2 ζ 2 2 , (28) where we introduced C = C γ γ . On the right hand side of the tangential force balance the active surface stress appears which we here consider symmetric and isotropic, i.e., ζ α β = ζ ( s 1 , s 2 ) δ α β . (29)

We note that the incorporation of anisotropic active stress is straightforward and illustrated in Ref. [42], which can easily be followed to include an anisotropic active stress evolving with the deforming cortex in the present algorithm. The derivatives of the active stress ∇ 1 ζ 1 1 = ∇ 1 ζ 2 2 = ∇ 1 ζ ( s 1 , s 2 ) , ∇ 2 ζ 1 1 = ∇ 2 ζ 2 2 = ∇ 2 ζ ( s 1 , s 2 ) , can be calculated from the given active surface stress distribution ζ(s ¹, s ²) on the deforming thin shell. From Eq. 29 we can directly use ∇ 1 ζ 1 2 = 0 , ∇ 1 ζ 2 1 = 0 ∇ 2 ζ 1 2 = 0 , ∇ 2 ζ 2 1 = 0 .

In our numerical implementation we discretize the thin shell representing the combination of cortex and membrane by a set of nodes, which we refer to by their position r _ν and/or index ν = 0, …, N−1, where N is the total number of nodes. The nodes are connected to flat triangles [65–68], as sketched in Figure 1. These triangles serve as a tool to determine the neighborhood of each node r _ν . For each node on the membrane we define a local coordinate system as in Ref. [42] and sketched on the right hand side of Figure 1. First, by averaging all normal vectors of the adjacent triangles weighted by angle [69], we obtain the local normal vector n on the node r _ν . By connecting the ν-th node to one of its neighbors and subtracting the normal component of the resulting vector, we determine the first in-plane coordinate vector e _ξ . The choice of the neighboring node is arbitrary as we show in Supplementary Appendix S3. By a normalized cross-product of n and e _ξ the second in-plane coordinate vector e _η is obtained. Thus, we have for each node a local coordinate system e ξ , e η , n ν = e ξ , e η , n . (30)

The local curvature tensor at the position of the central node expressed in the local coordinate system (Eq. 30) can be obtained by a parabolic fitting procedure with respect to node position as detailed in Ref. [42].

The active surface stress of the cortex is considered to be known, as detailed above, and its distribution can be described by an analytically tractable form on the initial shape of the cortex. As detailed in Ref. [42] the active surface stress distribution is therefore prescribed on the initial shape of the cortex for each node and afterwards advected with each node.

The key quantity of interest, the velocity v , is a three-dimensional vector in the three-dimensional Cartesian space v = v x e x + v y e y + v z e z , (31) with e _i , i = x, y, z the Cartesian unit coordinate vectors. Because the velocity field is defined on the thin shell, we can evaluate the velocity at each node ν giving v ν i = v i ( r ν ) . In total, the velocity field consists of the velocity of each node v x = { v ν x } ν ∈ { 0 ; … ; N − 1 } , v y = { v ν y } ν ∈ { 0 ; … ; N − 1 } , v z = { v ν z } ν ∈ { 0 ; … ; N − 1 } . (32)

The velocity components in the local coordinate system Eq. 30 can easily be obtained by projection onto the corresponding local coordinate vectors according to Eq. 6 v ξ = v ⋅ e ξ , v η = v ⋅ e η , v n = v ⋅ n . (33) With the derivatives of the velocity vector as introduced in Eq. 18 being evaluated in local coordinates, i.e., A _v = ∇ _ξ v , …, E _v = ∇ _ξ ∇ _η v , we can write down the velocity gradient Eq. 10 in the local coordinate system, e.g., the component v ξ η = 1 2 A v ⋅ e η + B v ⋅ e ξ .

For the actual calculation of the derivatives of the velocity vector we refer to section 3.5. Using the derivatives of the velocity and the velocity gradient in local coordinates, we are able to write down the force balance Eqs 16, 17 for each node ν in its local coordinate system, which corresponds to Eqs 26–28. Considering the tangential force balance Eqs 26, 27 together with the normal force balance Eq. 28, we end up with three differential equations. These we refer to in the following by its left hand side, which depends on the velocity vectors v _ν , and by the right hand side, which does not depend on the velocity vectors. Therefore we have three coupled differential equations which the velocity field has to fulfill l.h.s. 1 ν = r.h.s. 1 ν , l.h.s. 2 ν = r.h.s. 2 ν , l.h.s. n ν = r.h.s. n ν , (34) where the index 1,2 labels the two tangential force balance equations and the index n labels the normal force balance equation. The superscript ν emphasizes that we refer to the equation in local coordinates of node ν, which is the force balance evaluated at position r _ν of node ν.

Since only derivatives of the velocity field occur in the force balance equations, an arbitrary constant velocity can be added without violating the force balance equations. Furthermore, these derivatives are directional and under certain circumstances, e.g., axisymmetry, a velocity in azimuthal direction can be added. Therefore, we need to specify certain constraints to the velocity field, which correspond to boundary conditions for the partial differential equations.

First, we use a vanishing total normal velocity v ⁿ which corresponds to an incompressible liquid inside the cell, i.e., ∫ S v ⋅ n d S = ∑ ν = 0 N − 1 v ν ⋅ n A ν = 0 , (35) where the first equality results from the discretization of the integral into a sum over all nodes ν with A _ν being the local area of node ν. The local area is calculated using Meyer’s mixed area [42, 69]. This constraint of incompressibility accompanies the pressure difference P introduced in the normal force balance Eq. 17.

Second, we use that the total velocity of the cortex vanishes, i.e., ∫ S v d S = ∑ ν = 0 N − 1 v ν A ν = 0 . (36)

This constraint resolves the fact that due to the lack of an external force acting on the cortex the solution of the force balance equations is invariant under solid translation. In case of cell movement, one has to consider an additional coupling to the environment leading to additional external forces, for example via a friction coefficient in case of amoeboid motion [70] and in that case this constraint is lifted.

Third, we use with analogous motivation the fact that the effective total angular momentum of the cortex vanishes, i.e., ∫ S r × v d S = ∑ ν = 0 N − 1 r ν × v ν A ν = 0 , (37) in order to tackle the invariance under solid rotation. Again, in presence of additional external torques this constraint has to be lifted.

The overall goal is to solve the force balance equations for the velocity field. Accompanied by the incompressibility constraint from the cytoplasm, we also solve for the pressure. Since the solution fulfills the force balance on the complete thin shell, the force balance equations must be fulfilled at every node r _ν . Using Eq. 34, we define χ 2 = ∑ ν ( l.h.s. 1 ν − r.h.s. 1 ν ) 2 + ( l.h.s. 2 ν − r.h.s. 2 ν ) 2 + ( l.h.s. n ν − r.h.s. n ν ) 2 , (38) taking the two tangential and the normal force balance equation into account, with the sum ∑ ν over all nodes. The χ ² vanishes for the exact solution of the force balance equations by construction.

In order to incorporate the constraints detailed in the previous section 3.2, we use the concept of Lagrange multipliers. The incompressibility constraint Eq. 35 imposes one condition whereas the constraints in Eqs 36, 37 each impose three conditions due to the three components of the velocity vector. For each of the conditions we introduce a separate Lagrange multiplier. The constraint times the Lagrange multiplier is added to the χ ² in Eq. 38, e.g., χ 2 = ∑ ρ ∑ ν l.h.s. ρ ν − r.h.s. ρ ν 2 + λ 1 ∑ ν = 0 N − 1 v ν x A ν − 0 + λ 2 ∑ ν = 0 N − 1 v ν y A ν − 0 + λ 3 ∑ ν = 0 N − 1 v ν z A ν − 0 , (39) for the total velocity constraint Eq. 36 as illustration and with ρ = 1, 2, n.

For the exact solution the χ ² would become zero, numerically we seek the minimum of χ ². Thus, in simulations we have to solve for the velocity set v ν x , v ν y , v ν z in Eq. 32 together with the pressure difference P, which minimizes χ ². We therefore have to calculate arg min { v ν } , P χ 2 = arg min { v ν x } , { v ν y } , { v ν z } , P χ 2 . (40)

In order to perform this minimization, we calculate the derivatives of χ ² with respect to the velocity components of each node v ν i and with respect to P. Due to minimization the derivatives equal zero and we can separate this equality for terms linear in v ν i , P and terms independent of the velocity or pressure, where we use the parabolic expansion of the velocity field in Eq. 18. In total, this leads to a system of linear equations in matrix form with dimensions (3N + 1) × (3N + 1). There we have 3N rows due to the derivatives w.r.t. to the N velocities with 3 components each and we have 3N columns due to the N velocities with 3 components each in the solution vector. In addition, we have a row and a column due to the derivative with respect to and the dependency on the pressure, respectively. In the minimization procedure in Eq. 40, each Lagrange multiplier enters as an additional unknown, thus adding one row (containing derivatives w.r.t. the Lagrange mulitplier ∂ χ 2 ∂ λ j ) and one column containing the coefficients of λ _j .

In order to calculate the corresponding matrix to the system of linear equations we have to evaluate the left hand side and right hand side of all three force balance equations on each node in the local coordinate system Eq. 30. This we do using the velocity derivatives in local coordinates as detailed in 3.1. For details on the dependency on the actual velocity values we refer to section 3.6. The system of equations can finally be solved numerically using LU-Decomposition.

The velocity field in the cell cortex advects the cortex material and thus determines its reorganization and changes in shape over time. While tangential flows are parallel to the current shape, the normal velocity leads to a deformation of the cortex shape. For reasons of numerical stability, we consider here an additional damping contribution to the surface stress tensor which is similar to a bending viscosity [40]. The aim of this additional bending contribution, which is chosen to be small, is to include a condition on the second derivative of the normal velocity, which does not appear in the force balance equations. Due to its smallness and the limitations of the parabolic expansion we neglect the contribution of the bending contribution to the tangential force balance equations. The normal component of the velocity is obtained from the three-dimensional velocity by projection onto the normal vector Eq. 6. Derivatives are again obtained using a parabolic expansion similar to Eq. 18. We consider the time derivative of the curvature tensor in the Eulerian frame according to Salbreux and Jülicher [40]. ∂ t C α β = − ∇ α ∇ β v n − v n C α γ C β γ , (41) and the derivative of the trace ∂ t C γ γ = − ∇ γ ∇ γ v n − v n C γ δ C γ δ = − ∇ γ ∇ γ v n − v n C γ δ C δ γ . (42)

Using the expressions above, an additional in-plane surface stress related to a bending η ̄ and bending bulk η ̄ b viscosity, which account for dissipation due to curvature changes [40], is chosen as t b α β = 2 η ̄ g α γ ∂ t C γ β − 1 2 C δ δ δ γ β − C δ δ v n C γ β + η ̄ b ( ∂ t C γ γ ) g α β = 2 η ̄ − ∇ α ∇ β v n − v n C γ α C β γ − C γ γ v n C α β + η ̄ b − η ̄ − ∇ γ ∇ γ v n − v n C γ δ C δ γ g α β .

In Ref. [40] a contribution due to effective bending viscosities is introduced for a thin shell with broken up-down-symmetry and the bending viscosity can be negative. Again, we note that its addition here, in our three-dimensional framework, with an additional factor involving the normal velocity is motivated by reasons of numerical stability. The normal force balance contribution is calculated by contraction with the curvature tensor t b α β C α β = 2 η ̄ − ∇ α ∇ β v n − v n C α γ C β γ − C γ γ v n C α β C α β + η ̄ b − η ̄ − ∇ γ ∇ γ v n − v n C γ δ C δ γ g α β C α β = 2 η ̄ − ∇ α ∇ β v n − v n C α γ C β γ − C v n C α β C α β + η ̄ b − η ̄ − ∇ γ ∇ γ v n − v n C γ δ C δ γ C , with C = g ^αβ C _αβ . With the derivatives of the normal velocity A v n , B v n , … , E v n in analogy to Eq. 18, the trace of the second derivative can be written as ∇ γ ∇ γ v n = g α β ∇ α ∇ β v n = g 11 C v n + g 22 D v n . (43)

We further use C α β C α β = C 11 C 11 + 2 C 12 C 12 + C 22 C 22 = C ̃ , (44) and C α γ C β γ C α β = C ̄ , to obtain in total t b α β C α β = − 2 η ̄ C 11 − η ̄ b − η ̄ C g 11 C v n + − 2 η ̄ C 22 − η ̄ b − η ̄ C g 22 D v n − 4 η ̄ C 12 E v n + − 2 η ̄ C ̄ − 2 η ̄ C C ̃ − η ̄ b − η ̄ C C ̃ v n . (45)

This contribution can be added to the normal force balance equation and can be incorporated into the solution procedure detailed above in a straightforward manner.

Up to now, we have written the force balance Eqs 26–28 in terms of derivatives of the velocity vector in local coordinates A _v , B _v , …, E _v , which themselves are three-dimensional Cartesian vectors. What remains is to express these derivatives by means of the actual velocity vectors at all nodes, which allows us to write down the force balance as system of equations that is directly solvable for the velocity. To do so, i.e., to discretize the derivatives on the numerical grid, we use a parabolic expansion. For this, we consider the neighborhood of a node as sketched in Figure 1. Considering all N _ν neighboring nodes a(ν) of node ν, the coefficients A _v , B _v , …, E _v can be expressed by the velocity vector evaluated at the central node and at the sites of the neighboring nodes A v = A v v ν , { v a ( ν ) } with a ( ν ) = 1 , … , N ν . (46)

In order to do so we write down the velocity at some distance around the (central) node ν in the local coordinates of the central node ξ, η v ̄ ξ , η = v ν + A v ξ + B v η + 1 2 C v ξ 2 + 1 2 D v η 2 + E v ξ η , (47) which is the parabolic expansion Eq. 18 evaluated in the local coordinate system Eq. 30. The position vector can be evaluated at the position of the N _ν neighboring nodes in the local coordinates of the central node, i.e., ξ _a(ν), η _a(ν), which are obtained by projection of the difference vector between the two nodes onto the local coordinate vectors. The expanded velocity at the position of the neighbors then has to be equal to the actual velocity v _a of the neighbor node a(ν). Therefore, the squared difference of expanded and actual velocity has to be minimal χ v 2 = ∑ a = 1 N ν v ̄ ξ a ( ν ) , η a ( ν ) − v a 2 . (48)

The minimization is performed analytically in several steps, which are outlined here and detailed in the Supplementary Appendix S1. Using the differential form of Eq. 48 we calculate the derivatives with respect to A _v , B _v , …, E _v that have to equal zero due to minimization. This results in a system of linear equations expressed in terms of a 5 × 5 matrix for each of the three vector components with the velocity derivatives forming the solution vector. We invert the matrix analytically using Mathematica. By means of the inverse matrix we are able to write down the velocity derivatives A _v , B _v , …, E _v in terms of the velocity values, as illustrated by Eq. 46.

Putting everything of the above together, we here summarize our numerical algorithm. Starting from a prescribed active surface stress distribution ζ α β on the initial shape of the cortex, we compute the velocity in the cortex in Cartesian coordinates v by the following steps.

1) evaluate the active surface stress ζ α β on the nodes of the discretized cortex and construct a local coordinate system for each node.

As a first test setup, we compare three-dimensional numerical results on a static spherical cortex to corresponding analytical solutions. We consider an isotropic active surface stress distribution across the cortex ζ α β ( θ , ϕ ) = ζ ( θ , ϕ ) g α β (49) which is analytically expanded in terms of spherical harmonics Y _lm (θ, ϕ) ζ α β ( θ , ϕ ) = ζ ( θ , ϕ ) g α β = ∑ l = 0 ∞ ∑ m = − l l ζ l m Y l m ( θ , ϕ ) g α β , (50) where ζ _lm are the expansion coefficients of the active surface stress distribution. Here, greek indices denote the two angles of spherical coordinates θ and φ. For the given active surface stress profile in Eq. 50, we derive an analytical solution for the velocity field v in Supplementary Appendix S2. We expand the tangential velocity in terms of vector spherical harmonics v t = ∑ l , m v l m 1 ∂ α Y l m e α with coefficients v l m 1 (we note here that the component in direction n × ∂ α Y l m e α vanishes as discussed in the Supplementary Appendix S2) and the normal velocity in terms of spherical harmonics v n = ∑ l , m v l m n Y l m with coefficients v l m n . The force balance equations Eq. 16, Eq. 17 then lead to analytical relations for the velocity expansion coefficients v l m 1 and v l m n depending on ζ _lm . These relations are given in Supplementary Equations S2.23, S2.25 and determine the velocity field. In the following, we consider the quantities in units of the inital cortex radius R ₀ and the planar shear viscosity thus fixing R =R ₀=1 and η _s = 1. In the following, we consider the results for l >1. On the one hand, considering the full system of equations (case i)), the difference of the conditions Supplementary Equations (S2.23) and (S2.25) implies a vanishing in-plane velocity, i.e., the expansion coefficients of the in-plane velocity vanish v l m 1 , i = 0 . (51)

In turn, the normal velocity has to fulfill Supplementary Equation S2.27 which leads to v l m n , i = − R 2 η b ζ l m . (52)

In summary, the solution for the velocity field for given active surface stress ζ _lm is v i = − R 2 η b ζ l m Y l m n , (53) with vanishing tangential component.

On the other hand, this system can be utilized to test the in-plane velocity separately (case ii). To do so, we consider vanishing bulk viscosity η _b = 0 and fix the normal velocity v ⁿ to zero. The latter implies v l m n , i i = 0 , (54) instead of solving the normal force balance. In this case, the tangential force balance equations lead to the condition (Supplementary Equation S2.28) for the tangential velocity coefficients which become v l m 1 , i i = 1 l ( l + 1 ) − 2 ζ l m . (55)

Thus, for solving the tangential force balance equation only in case of a vanishing normal velocity, we obtain an analytical solution for the velocity v i i = v l m 1 , i i ∂ α Y l m e α = 1 l ( l + 1 ) − 2 ζ l m ∂ α Y l m e α , (56) which is due to a given active surface stress with coefficient ζ _lm .

In both cases, the analytical solution in Eq. 53 for 1) the full system and in Eq. 56 for 2) the pure tangential system allows a direct comparison between the analytical and numerical solution. Using the three-dimensional algorithm developed above, we consider a discrete cortex with 2562 nodes and 5120 triangles without bending viscosity. After applying an active surface stress distribution in terms of spherical harmonics according to Eq. 50, we solve numerically for the three-dimensional velocity field. In case 2) the normal force balance equation entering the minimization ansatz can be replaced by the condition of vanishing normal velocity on each node. In the following, the comparison to the analytical solution is performed by considering two types of error measures. On the one hand, we consider the absolute value of the difference between numerical v and analytical solution v ^an per node and average over all nodes ⟨ ϵ ⟩ = 1 N ∑ ν | v − v an | . (57)

On the other hand, we calculate the sum of the squared difference of the velocities relative to the maximum velocity, average, and take the square root to obtain the relative error ϵ rel = 1 N ∑ ν v x − v x an 2 + v y − v y an 2 + v z − v z an 2 | | v an, max | | . (58)

We first consider an axisymmetric active surface stress distribution in terms of ζ ₂₀ Y ₂₀ with ζ ₂₀ = 1 as shown in Figure 2A) on the discrete cortex. We further use η _b = 1 and an active surface stress offset of ζ ₀₀ = 0. A variation of the active surface stress offset leads to a finite pressure difference, but the velocity field is not altered. Figures 2B–E shows the full, three-dimensional velocity profile obtained numerically with the velocity magnitude per node color coded and the velocity direction indicated by arrows and the velocity depending on the polar angle θ in comparison to the analytical solution for both systems. In system 1), the normal velocity is finite and points outward the cell around the equator and inward at both poles. In system 2), the tangential velocity is finite and directed from the equator towards the poles. Each three-dimensional simulation shows very good agreement with the analytical solution, as shown by the comparison over the polar angle. As a next step, we quantify the error ⟨ϵ⟩ as well as ϵ _rel given in Eqs 57, 58, respectively, and vary the resolution of the discrete spherical cortex. As shown in Figures 2F,G, we obtain a systematic decrease of both errors with increasing resolution for both systems. The error shows a scaling inversely proportional to the number of nodes.

We then proceed in Figure 3 with a non-axisymmetric distribution of active surface stress in terms of ζ ₂₁ Y ₂₁ with ζ ₂₁ = 1 and ζ 00 = 4 π shown in Figure 3A. Again, we show the full system i) in Figures 3B,C and the purely tangential system ii) in Figures 3D,E. As a consequence of the non-axisymmetric active surface stress distribution, the cortical velocity profile is no longer axisymmetric. Figure 3B shows the three-dimensional velocity field for the full system i) and Figure 3D shows the three-dimensional velocity field for the purely tangential system ii) with normal velocity restricted to zero. Figure 3B shows four patches of large normal velocity. Opposite patches show either an outward or an inward pointing normal velocity. In-between, the normal velocity becomes zero. At the sites where in Figure 3B the normal velocity is maximal, the tangential velocity in Figure 3C vanishes. Similar to Figure 2 the tangential velocity points from sites with outward pointing normal velocity towards sites with inward pointing normal velocity in Figure 3B. In Figures 3C,E we compare the numerically obtained velocity to the analytical solution by showing the velocity value depending on the polar angle θ for an azimuthal angle ϕ = 0. Again, we observe very good agreement between numerical and analytical solution and a systematic convergence with increasing number of nodes as shown in Figures 3F,G for both systems i) and ii), respectively.

Having tested the algorithm for the static scenario in the previous section, we now apply our method to a dynamically evolving viscous active cortex immersed in a passive outer fluid. We consider an active surface stress distribution similar to what is known for cytokinesis, i.e., the separation of the two daughter cells by an evolving cleavage furrow during cell division [30, 32, 33, 46]. Initially, on the undeformed spherical cortex, we consider an active surface stress that is isotropic with constant offset ζ ₀ at the poles and increases towards the equator to form a contractile ring ζ 1 1 ( θ , ϕ ) = ζ 2 2 ( θ , ϕ ) = ζ ( θ , ϕ ) = ζ 0 + ζ ̂ exp − σ θ − π 2 p , (59) where ζ ̂ is the amplitude of active surface stress increase around the equator, p the exponent, and σ the width of the active surface stress distribution. Again we use R = 1 and for the active surface stress we fix ζ ₀ = 1, σ = 10, retaining the planar shear and bulk viscosity as η _s = η _b = 1. The exponent p is varied. We use an initially spherical mesh with 2562 nodes and 5120 triangles.

For the given active surface stress distribution in Eq. 59 and for given cortex shape, we solve numerically for the three-dimensional velocity field v on the discrete cortex. We split up the total velocity field v into the tangential velocity v t = 1 − n ⊗ n ⋅ v and the normal velocity v n = n ⊗ n ⋅ v . Since only the normal velocity leads to flows in the cortex that change the overall cortex shape, we integrate the cortex evolution similar to Refs. [45, 71] by advecting each node ν with its normal velocity v n ν . With this, a numerically demanding re-meshing, which would be necessary in case of tangential advection of the mesh nodes, is avoided. Advection is carried out using the Euler algorithm r ν ( t + Δ t ) = r ν ( t ) + v n ν ( t ) Δ t , (60) with a time step Δt, which is typically chosen to be Δt = 5 × 10^–4 t _a with the time scale t a = η s ζ 0 . After a number of time steps, we observe that the normal velocity on the whole cortex vanishes, i.e., the cortex has reached a final, steady state.

We first consider a broad distribution of active surface stress with exponent p = 8 in Eq. 59 with ζ ̂ = 1 as shown in Figure 4A. We apply a bending viscosity of η ̄ = 0.0005 η s and η ̄ b = 0.0025 η s . Figure 4B shows the velocity field on the initially spherical cortex. Because of the dominating contractile active surface stress around the equator, the normal velocity points towards the cell interior around the equator. As a consequence of the incompressibility of the enclosed fluid, the normal velocity points outward at the poles. Therefore, the velocity leads to an expansion of the cortex along the z-axis and a simultaneous contraction around the equator. Integration of the cortex evolution leads to the shape shown in Figure 4C, where we also show the final velocity field on the deformed cortex. Here, no normal velocity remains, but a finite tangential velocity points towards the equator.

For this setup, we compare the three-dimensional dynamics to the one obtained by axisymmetric simulations of viscous active surfaces [72] in Figure 5. Details of the axisymmetric simulations are given in Supplementary Appendix S4. For comparison between the simulation methods, we track the distance between the poles, which is divided by a factor of two, as well as the radius of the furrow at the equator. Both quantities are displayed as functions of time in Figure 5A. While the pole to pole radius increases over time, the cortex contracts at the equator and thus the furrow radius decreases. Both quantities reach a constant plateau at longer time, which corresponds to the convergence of the simulation and results from the interplay of the contractile active stress at the pole and cortical ring contraction mediated by the incompressibility of the enclosed fluid [33]. We obtain excellent agreement between three-dimensional and axisymmetric simulation.

To go one step further, we also compare the velocity field on the discrete cortex at different times. On the left hand side of Figure 5 we show the absolute value of the tangential velocity v _t depending on the axial position z. On the right hand side we show the normal velocity v _n . Velocities are shown at different times and are normalized by v a = R 0 t a . While the tangential velocity increases over time, i.e., with increasing deformation of the cortex, the normal velocity decreases until the cortex reaches its final shape (the latter is shown in Figure 4C). Overall, our developed three-dimensional algorithm leads to a velocity field on the evolving cortex, which is in very good agreement with the axisymmetric simulation.

Next, we consider an active surface stress distribution around the equator, which is narrower with exponent p = 4 in Eq. 59 and systematically vary the amplitude ζ ̂ . Figure 6 shows the active surface stress distribution for ζ ̂ = 1 as illustration in Figure 6A together with the corresponding shape and velocity profile in the Figure 6B initial and Figure 6C final state. Figures 6D,E show the initial and final state for ζ ̂ = 3 . We use η ̄ = 0.0005 η s and η ̄ b = 0.0025 η s in Figures 6A–C and η ̄ = 0.0001 η s , η ̄ b = 0.0035 η s in Figures 6D,E. With increasing amplitude of the active surface stress the cortex extends farther at the poles and constricts around the equator more strongly. The final tangential velocity increases systematically, which is due to the gradient in the active surface stress distribution increasing with increasing magnitude and constant offset. We again compare the dynamics based on the pole to pole and the furrow radius as well as the velocity field at different times to axisymmetric simulations. In Figure 7, we show the comparison for the amplitude ζ ̂ = 1 and in Figure 8 for ζ ̂ = 3 . In all cases the results of the three-dimensional algorithm agree very well with the axisymmetric simulations.

Overall, the comparison to axisymmetric simulations clearly shows the accuracy of the developed algorithm for a viscous active cortex in three dimensions also in case of dynamic cortex evolution.

As a first sample application we consider a non-axisymmetric stress pattern, which leads to the formation of a local indentation reminiscent of fold formation during embryogenesis [50]. We apply a line of increased active stress connecting both poles, which broadens towards the equator. Therefore, we modulate the active stress pattern in Eq. 59 with two additional Gaussian functions in ϕ-direction, i.e., we multiply by exp − σ ϕ ( ϕ − π / 2 ) 4 or exp − σ ϕ ( ϕ − 3 π / 2 ) 4 , respectively. The resulting active surface stress distribution is shown in Figure 9A for ζ ̂ = 3.5 and σ _ϕ = 200 and in Figure 9D for ζ ̂ = 4.5 and σ _ϕ = 12, where in both cases ζ ₀ = 1, p = 4, and σ = 1.5. We note that the cortex fulfills up-down symmetry such that we only show one half for clarity. The velocity profiles shown in Figures 9B,C,E,F, respectively, show a large initial inward velocity across the patch of enhanced active stress, where in the end Figures 9C,F a local fold forms. During fold formation an in-plane flow from the poles towards the fold develops. At the poles the cortex expands.

In order to provide a further example application of an evolving cortex in a non-axisymmetric situation, we consider in the following a cell cortex with an initial shape that is not aligned with respect to the Cartesian coordinate axes. We use the same values for the viscosity as above, i.e., η _s = 1 and η _b = 1, an exponent p = 4 for the active surface stress distribution with constant offset ζ ₀ = 1 and amplitude ζ ̂ = 3.5 . For the bending viscosity we choose η ̄ = 1 × 1 0 − 4 and η ̄ b = 5 × 1 0 − 3 . Initially, we deform the cortex according to a shear strain of γ _x = 0.35 along the x-axis.

For the first setup, the active surface stress distribution is fully non-axisymmetric as shown in Figure 10A. We show in Figures 10B–E the dynamic evolution of the cortex with the velocity magnitude per node color coded and the velocity direction indicated by arrows. Initially, a large normal velocity is observed in Figure 10B with a stripe-like pattern. Roughly along the long axis of the sheared cortex the velocity points inwards the cell cortex. In direction perpendicular to the shear, the velocity points outwards of the cell cortex. This illustrates the tendency of the cell cortex to realign with respect to the symmetry of the active surface stress distribution on the undeformed cortex, which is similar to the one shown in Figure 6A. The reorientation is visible throughout the time evolution as well. Finally, the cell cortex is nearly aligned with the initial gradient of the active surface stress distribution. This points to the symmetry of the active surface stress distribution being the driving mechanism behind this reorientation. From Figures 10D,E it can be seen that the cortex subsequently contracts around the equator at the end of the reorientation. In total, this indicates that the underlying symmetry of the active surface stress distribution triggers a reorientation of the total cortex, which in the end leads to a robust furrow formation despite the initial shear disturbance. We hypothesize that in nature this autonomous reorientation may explain the robustness of cell division in noisy environments where external non-axisymmetric perturbations due to, e.g., forces by neighboring cells may be prevalent.