A Three-Dimensional Numerical Model of an Active Cell Cortex in the Viscous Limit

The cell cortex is a highly dynamic network of cytoskeletal filaments in which motor proteins induce active cortical stresses which in turn drive dynamic cellular processes such as cell motility, furrow formation or cytokinesis during cell division. Here, we develop a three-dimensional computational model of a cell cortex in the viscous limit including active cortical flows. Combining active gel and thin shell theory, we base our computational tool directly on the force balance equations for the velocity field on a discretized and arbitrarily deforming cortex. Since our method is based on the general force balance equations, it can easily be extended to more complex biological dependencies in terms of the constitutive laws or a dynamic coupling to a suspending fluid. We validate our algorithm by investigating the formation of a cleavage furrow on a biological cell immersed in a passive outer fluid, where we successfully compare our results to axi-symmetric simulations. We then apply our fully three-dimensional algorithm to fold formation and to study furrow formation under the influence of non-axisymmetric disturbances such as external shear. We report a reorientation mechanism by which the cell autonomously realigns its axis perpendicular to the furrow plane thus contributing to the robustness of cell division under realistic environmental conditions.

Materials containing cytoskeletal filaments and motor proteins are successfully described by active gel theory [2,[35][36][37]. Its key ingredient is the actively induced force by the motor proteins, which leads to an active stress in the material [38,39]. In addition, polymerization and depolymerization lead to transient changes within the active gel [20,35]. Focusing on the cell cortex of actomyosin or on epithelial tissue, active gel theory has recently been formulated in the framework of thin shell theory [40]. The generic active gel theory [35] incorporates the viscoelastic nature of the cytoskeletal assemblies, where the presence of an active stress triggers both elastic deformations and viscous flows. Starting from the viscoelastic theory two limits of the cortex behavior can be considered, the elastic limit [40][41][42][43] and the viscous limit [21,28,30,33,40,44,45]. In the latter, viscous flows arise which are responsible for a large number of cellular processes. In cell morphogenesis transitions from a spherical initial shape to pear-like [20,21] or oblate [8] shapes have been reported. Cell division has successfully been modeled [30,31,33,34,44,46] including a threshold active stress needed to complete cytokinesis [33]. Further examples include cell movement [24,47] and embryogenesis [19,48,49].
Such investigations involving the cell cortex are based on either analytical calculations [21,30,41,44] or numerical simulations [9,33,42,45,47,50,51]. Farutin et al. [47] investigated cell crawling by coupling cortical mechanics to the boundary integral equation in an axisymmetric treatment. Turlier et al. [33] used a numerical, axisymmetric formulation of cell cytokinesis advecting tracer points to track the position of the membrane. Mietke et al. [9] developed another axisymmetric simulation method by discretizing the arc length of the cellular membrane. Including myosin activity in terms of a preferred curvature in the framework of an elastic shell, Heer et al. [50] determined the equilibrium of a tissue shell in three dimensions. Torres-Sánchez et al. [45] developed a fully three-dimensional computational model of a viscous active cell cortex in the absence of an external fluid. Being based on Onsager's variational principle [52] their method utilizes a subdivision finite element scheme for discretization. A generalization of the latter to arbitrary topology using local Monge parametrization method has been provided [51]. Recently, two of us have developed a dynamic three-dimensional simulation model of an active cell cortex in the elastic limit based on a parabolic fitting procedure [42]. This approach provides the advantage that it is dynamically coupled to the surrounding fluid flow using the lattice-Boltzmann/immersed boundary method and thus allows one to study cortical dynamics in realistic environments such as, e.g., blood flow [43].
In the present manuscript we develop a new method, transferring the ideas of ref. [42] to the viscous limit in order to numerically obtain the velocity profile in a three-dimensional, thin, and arbitrarily deformed active cell cortex. We use the thin shell formulation combined with active gel theory [40] to obtain the force balance equations for the cortex involving active and viscous stresses. In contrast to the variational approach of Refs. [45,51] we directly start with the force balance equations of the cortex, which are discretized using a parabolic fitting procedure. We furthermore base our approach on the velocity vector expressed in three-dimensional Cartesian coordinates. Using an analytical inversion of the parabolic expansion and fitting of the cortical velocity field, we evaluate the force balance equation on the discrete nodes of the cortex. Solving the resulting system of coupled equations globally on the cortex by means of a minimization ansatz, we finally solve for the cortical velocity field. Considering the normal component of the velocity, which characterizes the flow of the cortical material leading to changes of the cortex shape, the cortex shape can be evolved in time in order to obtain the deforming cell shape. In turn, the velocity field is obtained on the continuously evolving cortex. We provide an in-depth validation using analytical results on an undeformed, spherical cortex and axisymmetric simulations for the evolving shape and cortical flow field. As an application, we consider an initial shear deformation of the cortex and investigate its dynamic evolution. We find that the cell realigns its axis perpendicular to the furrow plane thus autonomously rectifying non-axisymmetric external perturbations. Due to its simplicity, our proposed algorithm can be the basis for a dynamic coupling of a viscous active cortex or a tissue to a suspending fluid using, e.g., a coupled immersed boundary and lattice-Boltzmann method.
We first introduce the force balance equations of the cell cortex treated as a viscous thin shell in section 2. Afterwards, we present the discretization of the cortex based on a parabolic fitting procedure and the numerical solution procedure in section 3. In section 4, we validate our method to analytical calculations and axisymmetric simulations in a detailed manner and apply our algorithm to dynamic fold and furrow formation in response to a non-axisymmetric shear deformation. Finally, we conclude in section 5.

THIN SHELL FORMULATION OF A VISCOUS ACTIVE CORTEX
In the following we consider a cell cortex in the viscous limit subject to internal active stresses, which lead to effective flows within the actomyosin network [37,40]. Since the cortex is typically very small compared to the cell diameter [1], it is considered as a thin shell [40,53], i.e., as a two dimensional manifold in three-dimensional space. The framework for the mathematical description of a thin shell is differential geometry [54,55] which we introduce in section 2.1. In section 2.2 we provide the analytical formulation of the force balance for the viscous active cortex, which we express in terms of the velocity and its derivatives in section 2.3.

Differential Geometry
In general, the two dimensional thin shell representing the cell cortex is parametrized by the three-dimensional vector X (s 1 , s 2 ) which depends on the two surface coordinates s 1 , s 2 . 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 and using those a unit normal vector on the thin shell that points outwards can be defined The metric tensor of the thin shell is defined by and its curvature tensor by with Greek indices α, β, c 1, 2 referring to the in-plane coordinates s 1 , s 2 . 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 αc g cβ . Correspondingly, an upper (contravariant) index can be lowered by a β g βα a α or t αβ t α c g cβ . 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 z α f and the covariant derivative by ∇ α f. The covariant derivative is defined for a scalar by ∇ α f z α f, for a vector component by ∇ α a β z α a β + Γ β αc a c , and for a general tensor by ∇ α t βc z α t βc + Γ β δα t δc + Γ c αδ t βδ , where Γ c αβ 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 βc ∇ c . The covariant derivative of the in-plane and normal coordinate vector is given by the equations of Gauss and Weingarten 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 1 , s 2 ), 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 α and v n v · n, with its in-plane components v α and its normal component v n such that v v α e α + v n n.
The total in-plane vectorial component of the velocity field can further be determined by 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].
which is equivalent to the definition in vector notation

Force Balance for a Viscous Active Cortex
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 inplane surface stress t αβ , a 2 × 2-tensor, and the normal surface stress t α n [42] as t α t αβ e β + t α n n.
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.
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].
with the viscous surface stress tensor determined by 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 −4 N/m [57][58][59][60][61]. A typical value of the shear viscosity of the cell cortex is η s 27 × 10 −4 Pa s m [62]. For the plasma membrane of a red blood cell a shear viscosity of about η s ≈ 3 × 10 −7 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.
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 Force Balance in Terms of Cartesian Vectors
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 1 , s 2 ) in the surrounding of a given position r ] on the thin shell v s 1 , 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 inplane 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 The trace of the velocity gradient is v c c v αβ g αβ g 11 v 11 + 2g 12 v 12 + g 22 v 22 Next, we can calculate the derivative of the velocity gradient tensor in Eq. 10 which, using the equations of Weingarten and Gauss (5), becomes Using that ∇ α g βc 0, we can calculate the derivative of the contravariant components of the velocity gradient by Furthermore, we calculate the derivative of the velocity gradient's trace via using the formula for the gradient in Eq. 22. Writing down these equations for fixed indices α, β, c 1, 2, they can be explicitly written in terms of 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 βc g cα and g αβ g βc δ c α leads to 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 the second tangential force balance equation and the normal force balance equation where we introduced C C 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., Frontiers in Physics | www.frontiersin.org October 2021 | Volume 9 | Article 753230 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 can be calculated from the given active surface stress distribution ζ(s 1 , s 2 ) on the deforming thin shell. From Eq. 29 we can directly use

COMPUTATIONAL METHOD FOR A VISCOUS ACTIVE CORTEX
As introduced above, our aim is a computational method to calculate the velocity field corresponding to a given active surface stress on a deforming cortex. In the following, we illustrate the steps of this numerical calculation. First, we discretize the thin shell representing a cell cortex and the velocity field in section 3.1. This allows us to write down the force balance equations Eq. 26, Eq. 27, Eq. 28 for each node of the discrete thin shell. We introduce constraints on the velocity field such as vanishing total velocity or angular momentum in section 3.2. In order to obtain the solution of velocity field we use a minimization ansatz as detailed in section 3.3. For a better numerical stability, we introduce an addition to the surface stress motivated by a bending viscosity in section 3.4. In order to perform the minimization ansatz and in turn to solve the system of equations we need the derivatives of the velocity field in Eq. 18 as functions of the velocity values at all discrete nodes as illustrated in section 3.5. The final step consists of building a matrix for the equation system as detailed in section 3.6.

Discretization of the Cortex
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][66][67][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 .
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 threedimensional vector in the three-dimensional Cartesian space 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 FIGURE 1 | Cortex discretization. Discretization of a thin shell by nodes connected to triangles (left). For each node r ν , a local coordinate system with in-plane coordinate vectors e ξ , e η and normal vector n is constructed. The velocity field, its derivatives and the force balance equations can be expressed in this local coordinate system (e ξ , e η , n) at the position r ν .
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 With the derivatives of the velocity vector as introduced in Eq. 18 being evaluated in local coordinates, i.e., we can write down the velocity gradient Eq. 10 in the local coordinate system, e.g., the component 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 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 ].

Constraints
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 n which corresponds to an incompressible liquid inside the cell, i.e., 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., 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., in order to tackle the invariance under solid rotation. Again, in presence of additional external torques this constraint has to be lifted.

Minimization Ansatz for the Force Balance
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 taking the two tangential and the normal force balance equation into account, with the sum ] over all nodes. The χ 2 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 χ 2 in Eq. 38, e.g., In order to perform this minimization, we calculate the derivatives of χ 2 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 zχ 2 zλ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.

Improving Numerical Stability by an Artificial Bending Viscosity
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 threedimensional 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].
and the derivative of the trace 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 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 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 We further use and to obtain in total 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.

Velocity Derivatives as a Function of Neighboring Velocity Vectors
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 , . . ., Frontiers in Physics | www.frontiersin.org October 2021 | Volume 9 | Article 753230 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.
In order to do so we write down the velocity at some distance around the (central) node ] in the local coordinates of the central 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 The minimization is performed analytically in several steps, which are outlined here and detailed in the Supplementary The three-dimensional velocity field is shown as obtained by the numerical solution of the full system. (C) While the tangential velocity is zero, the normal velocity depending on the polar angle θ agrees very well with the analytical prediction. (D) The three-dimensional velocity field is shown as obtained by the numerical solution for solving the tangential force balance only. (E) While the normal velocity is zero, the absolute value of the tangential velocity depending on the polar angle θ agrees very well with the analytical prediction. For both the full system i) in (F) and the purely tangential system ii) in (G), the error of the velocity converges with increasing resolution.
Frontiers in Physics | www.frontiersin.org October 2021 | Volume 9 | Article 753230 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.

Summary of the Algorithm
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. 2) compute derivatives of the active surface stress ζ β α together with the metric and curvature tensor at each node using a parabolic fitting for the active surface stress or the position vector, respectively, analogous to Eq. 18.

3) use the inverse parabolic fitting detailed in Supplementary
Appendix S1 to obtain the derivatives of the velocity vector up to second order A v , B v , . . ., E v in local coordinates at each

VALIDATION AND APPLICATION
In the following, we intensively validate the developed algorithm for a viscous active cell cortex in three dimensions. First, we consider in section 4.1 a static spherical cortex subjected to an active surface stress distribution in terms of spherical harmonics. We compare the resulting velocity field on the cortex to analytical solutions derived in Supplementary Appendix S2. Next, we compare our three-dimensional simulations for a dynamically deforming cortex to axisymmetric simulations in section 4.2. Finally, we exploit the capabilities of our fully 3D algorithm to study fold formation in section 4.3 as well as furrow formation under nonaxisymmetric external perturbations in section 4.4.

Velocity Profile by Spherical Harmonics on a Static Spherical Cortex
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 (θ, ϕ) 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 1 lm z α Y lm e α with coefficients v 1 lm (we note here that the component in direction n × (z α Y lm )e α vanishes as discussed in the Supplementary Appendix S2) and the normal velocity in terms of spherical harmonics v n l,m v n lm Y lm with coefficients v n lm .
The force balance equations Eq. 16, Eq. 17 then lead to analytical relations for the velocity expansion coefficients v 1 lm and v n lm 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 0 and the planar shear viscosity thus fixing R R 0 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 In turn, the normal velocity has to fulfill Supplementary Equation S2.27 which leads to v n,i) In summary, the solution for the velocity field for given active surface stress ζ lm is with vanishing tangential component.
On the other hand, this system can be utilized to test the inplane velocity separately (case ii). To do so, we consider vanishing bulk viscosity η b 0 and fix the normal velocity v n to zero. The latter implies 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 Thus, for solving the tangential force balance equation only in case of a vanishing normal velocity, we obtain an analytical solution for the velocity 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 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 We first consider an axisymmetric active surface stress distribution in terms of ζ 20 Y 20 with ζ 20 1 as shown in Figure 2A) on the discrete cortex. We further use η b 1 and an active surface stress offset of ζ 00 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, threedimensional 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 ζ 21 Y 21 with ζ 21 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 nonaxisymmetric 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.

Dynamics of an Axisymmetric Viscous Active Cortex
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 ζ 0 at the poles and increases towards the equator to form a contractile ring 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 ζ 0 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 threedimensional 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 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  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.

Embryonic Fold Formation
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

Robust Cleavage Furrow Formation Under Non-axisymmetric Perturbations
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 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 nonaxisymmetric perturbations due to, e.g., forces by neighboring cells may be prevalent.
In the second setup, we consider the dynamic evolution of the cortex subject to the same initial shear, but with an active surface stress distribution that is rotated such that the gradient of the active surface stress distribution nearly aligns with the shear. The initial shape is shown in Figure 11A while the evolution of the cortex is shown in Figures 11B-E at the same time points as above. Initially, the flow field shows four patches of large velocity in Figure 11B. The evolving flow field first leads to a broadening of the cortex and a slight reorientation, presumably with respect to the symmetry of the active surface stress distribution. Subsequently, the cortex contracts around the equator with the eventual shape shown in Figure 11E. The tangential velocity is directed from the poles to the equator. We note that for very long times-after reaching a steady state regarding the reorientation and furrow formation-a slight asymmetry in the volume of the two daughter cells leads to an unstable shape of the dividing cell, which is due to the difference in the Laplace pressure following a difference in curvature between both regions [31,46].
Finally, we consider about twice the shear rate, i.e., c x 0.75, in Figure 12 where the results are similar to before.
These two examples illustrate the applicability of our algorithm to non-axisymmetric situations with a threedimensional deformation of the cortex.

CONCLUSION
In this work, we presented a simulation algorithm for the timedependent and fully three-dimensional deformation of a viscous cell cortex including active stresses. The cell cortex is represented as a thin shell in the framework of active gel theory and discretized using a set of nodes connected by flat triangles. Using an inverted parabolic fitting procedure, we were able to express the governing force balance equations on each node in terms of the unknown velocity components. The resulting linear system is solved on the whole cortex using a global minimization ansatz. Extensive comparison to analytical solutions on a rigid sphere as well as to axisymmetric simulations of a dynamically deforming shape showed very good agreement. As first examples of non-axisymmetric situations, we considered fold formation as well as a cortex subject to a shear deformation showing a reorientation mechanism of the cortex, which can further contribute to robustness of furrow formation or cell division.
The inherent three-dimensional character of the present algorithm allows the computational investigation of more complicated scenarios of cell or tissue mechanics in the future. Especially, the robustness of cell division/furrow formation with respect to externally induced deformations, external constraints as e.g. present during embryogenesis, or even positioning of the contractile ring can be investigated in more detail. Shape changes of tissue layers triggered by active stress distributions as they occur during embryonic or cancer development can also be addressed.
The generality of the presented algorithm, which directly solves the force balance equations, allows for a straightforward inclusion of different constitutive laws for cortex behavior and active stress. By way of example, the active surface stress distribution on the discrete cortex could be derived from a concentration field and/or take into account cytoskeleton polymerization and depolymerization. Coupling to an external environment can be achieved by adding corresponding external forces to the force balance equations. Therefore, it would be possible to compute cortex dynamics inside a flowing liquid in combination with a fluid solver. Our developed numerical tool can thus be the basis for future investigations of various three-dimensional scenarios of cell and tissue morphogenesis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. An implementation of the developed three-dimensional algorithm is available on request to the corresponding author(s). The source code for the axisymmetric simulations that has been used for validation can be found online: https://github.com/DianaKhoromskaia/ AxisymCortex.

AUTHOR CONTRIBUTIONS
CB designed the research, performed, analyzed, and interpreted the results, and wrote the article. DK performed the research and wrote the article. GS designed the research and wrote the article. SG designed the research, interpreted the results, and wrote the article.

FUNDING
DK and GS were supported by the Francis Crick Institute which receives its core funding from Cancer Research UK (FC001317), the UK Medical Research Council (FC001317), and the Wellcome Trust (FC001317), and by a CRUK multidisciplinary project award (C55977/A23342). This publication was funded by the University of Bayreuth Open Access Publishing Fund. We acknowledge funding from Deutsche Forschungsgemeinschaft in the framework of FOR 2688 "Instabilities, Bifurcations and Migration in Pulsating Flow," Project B3 (417989940). CB thanks the Studienstiftung des Deutschen Volkes for financial support.