A viscous active shell theory of the cell cortex

The cell cortex is a thin layer beneath the plasma membrane that gives animal cells mechanical resistance and drives most of their shape changes, from migration, division to multicellular morphogenesis. It is mainly composed of actin filaments, actin binding proteins, and myosin molecular motors. Constantly stirred by myosin motors and under fast renewal, this material may be well described by viscous and contractile active-gel theories. Here, we assume that the cortex is a thin viscous shell with non-negligible curvature and use asymptotic expansions to find the leading-order equations describing its shape dynamics, starting from constitutive equations for an incompressible viscous active gel. Reducing the three-dimensional equations leads to a Koiter-like shell theory, where both resistance to stretching and bending rates are present. Constitutive equations are completed by a kinematical equation describing the evolution of the cortex thickness with turnover. We show that tension and moment resultants depend not only on the shell deformation rate and motor activity but also on the active turnover of the material, which may also exert either contractile or extensile stress. Using the finite-element method, we implement our theory numerically to study two biological examples of drastic cell shape changes: osmotic shocks and cell division. Our work provides a numerical implementation of thin active viscous layers and a generic theoretical framework to develop shell theories for slender active biological structures.


Introduction
Biological cells can exhibit a large diversity of shapes, which are intimately linked to the multiple key functions they perform: division, polarisation, apoptosis (i.e., programmed death), or migration. Dynamic changes in the shape of animal cells, also called morphogenesis, is mainly controlled by the cell cortex, a dense meshwork of filaments underlying the plasma membrane, capable of generating internal active tensions and torques [1].
The cell cortex is a thin biological structure of few hundreds of nanometers in thickness that extends over several tens of micrometers, the typical cell size [2]. It is composed of semi-flexible actin polymers, actinbinding proteins, and myosin molecular motors, which are the substrates for various phosphorylation events keeping this material out of equilibrium [3]. By continuously hydrolyzing ATP (adenosine triphosphate) into ADP (adenosine diphosphate) and inorganic phosphate, actomyosin networks have two remarkable properties, which are relatively unusual in classical mechanics. First, they are under permanent renewal, principally through the assembly and disassembly of filaments, which happens within seconds; this fast material turnover releases elastic stresses accumulated in the network, making this material viscous on timescales longer than a few dozens of seconds [4]. Second, they are the subject of internal active loads created by molecular motors, which exert contractile traction forces on neighboring filaments; this creates a source of active stress that is continuously dissipated and renewed, leading overall to an effective surface tension at the cell level [5,6,7,8]; it can also lead to structural rearrangements of actin filaments within the plane of the cell's surface, which has motivated the recent emergence of active-gel theories [9]. Relying on mass and momentum balance and exploiting the fundamental symmetries of the system, active-gel theories extend hydrodynamic theories for nematic or polar liquid crystals [10] to account for the activity of motors and, to a certain extent, for the renewal of the material [11,9]. They have provided a convenient generic theoretical framework to describe the large-scale mechanical properties of active materials through hydrodynamic variables. The cortex mechanics is consecutively described in general by two sources of stress: an active stress, which is controlled spatiotemporally by the cell metabolism and depends on the local filament organization, and a passive stress, which results from the viscoelastic deformation of the network in response to external loads and internal active stress. This theoretical framework has proven its ability to capture the essential physics characterizing the morphological changes observed experimentally in single cells [5,6,7,12,13,14,15,16].
Because of the slender character of the cell cortex (around 0.2-2 in thickness [17,2] for a typical cell radius of , it is reasonable to adopt thin shell theories to model its mechanics at the cell level [18,19]. Shells are considered thin if their thickness is much smaller -typically by one order of magnitude or more -than other characteristic dimensions, such as the curvature radius. The main goal of thin shell theories is to reduce the three-dimensional equations into an equivalent two-dimensional model, which is generally defined onto the midsurface (i.e. the surface that splits the shell into two layers of equal thickness at any point), although different choices of reference surfaces are possible [20,21]. To perform this dimensional reduction, two generic approaches are usually undertaken: one may calculate the limit, as the slenderness parameter goes to zero, of the the nonlinear 3D model using asymptotic methods; alternatively one may restrict, from the beginning, the range of admissible deformations and stresses used in the nonlinear 3D shell model, by means of -generally ad hoc -assumptions, for instance the Kirchhoff-love hypothesis [22]. Even though the passive response of the cortex is essentially viscous at the timescale of interest, we expect here our viscous active shell model to be closely related elastic plates and shell theories, which have been extensively studied in solid mechanics. This is due to the Stokes-Rayleigh analogy, which states that the equations of slow viscous flow and linear elasticity are interchangeable in the incompressible limit if velocities replace displacements and strain rates replace strains [23,24,25,26]. For a more detailed introduction on elastic shell theories, we refer to classic textbooks [27,28,29,30]. For viscous materials, only a limited number of thin sheet models have been previously derived. A few membrane models (in which bending resistance is neglected) have been developed in arbitrary shapes [31,32] but models accounting for bending resistance were mostly derived for specific geometries [33,34,35], with one exception in general curvilinear coordinates [36]. Relying on an asymptotic approach, this last study constitutes the starting point of our active shell theory for the cell cortex.
Most of the recent works on dynamic cell shape changes have considered the mechanics of the cortex, modeled as an active and viscous surface equipped with tensions only [37,7,38,39]. Other works have focused on the lipid bilayer (or plasma membrane) mechanics (onto which the cortex is attached through specialized linker proteins), which is generally described as a 2-dimensional fluid layer equipped with an elastic bending rigidity [40,41,42,43,44,45]. Lipid bilayers may indeed be regarded as a stack of two lipid monolayers of nanometric thickness that do not flow transversely, which justifies why viscous dissipation is inherently 2-dimensional in that case. Here, we argue that the lipid bilayer contribution to cell mechanics is generally negligible, and we focus on the cortex. The cortex present several features that justify the need to derive a proper thin viscous shell theory starting from 3-dimensional constitutive equations. First, the cortex is a bulk continuous material, with a thickness that may vary along the cell's surface, which is in contrast to the classical assumption of constant thickness in thin elastic shell theories [46] or lipid bilayers. Second, the turnover of cortical material is not homogeneous along the transverse sheet direction: polymerization happens mainly at the lipid bilayer surface, through membrane-attached actin nucleator proteins called formins [47], while depolymerization acts in bulk, idealistically at each filament ends [48]. As commonly described, we assume that this mode of turnover breaks the symmetry of the shell with respect to its midsurface by creating an inward flux of actomyosin toward the cell interior, as depicted schematically on Fig. 1. This flux can potentially generate tangential viscous flows of cortical material, and therefore and additional active source of stress, that has been neglected so far. Third, pure membrane models of the cell surface cannot handle compressive forces nor active torque generation, that may lead to surface buckling. Finally, if active-gel theories provide relevant bulk constitutive equations for cortical material, the mathematical form of equivalent two-dimensional stress and torque resultants after dimensional reduction is hard to guess readily. Figure 1: Schematic representation of the cell cortex, composed of actin filaments and actin-binding proteins, in particular myosin molecular motors. The cortex forms a thin layer of thickness ≪ 0 beneath the plasma membrane, where 0 is the typical curvature radius of the cell. It is under permanent turnover, with a polymerization flux from the plasma membrane surface and a depolymerization rate in the bulk of the layer. The dimensional reduction into a thin theory theory consists in finding the two-dimensional approximate in-plane and transverse shear stress and , and torque resultants, defined on the middle-surface of the shell.
This paper proposes therefore a general thin active viscous shell theory that aims to predict how the shape of a cell responds to arbitrary loads. The load can be internal, such as gradients of active stress, or external, originating from the mechanical interaction of the cell with its environment. Extending to active-gel materials a previous covariant description of thin viscous shells [36], we present a systematic asymptotic expansion of governing equations in powers of a slenderness parameter ≡ 0 ∕ , where 0 and are respectively the characteristic thickness and a length scale over which the flows varies across the sheet (for a cell, is typically its radius 0 , and the thickness is fixed by turnover 0 ∼ , see Fig. 1). We start from three-dimensional active-gel constitutive laws and mass balance equations. The expansion is performed in two limits, stretching and bending-dominated, which are then combined to build an approximate two-dimensional model that is valid for the whole range of intermediate behaviors, following the generic method of matched asymptotic expansion [49]. The two-dimensional model is complemented by a kinematic equation describing the evolution of the thickness. Aiming at facilitating the spread, modification, refinement of our model and to highlight its features, we furthermore devise a numerical scheme [50] (available on https://github.com/VirtualEmbryo/active_viscous_shell) to solve the 2-dimensional thin shell equations with the finite-element library FEniCS [51], and we illustrate numerically two drastic dynamic cell shape changes: cell division and osmotic shocks.
The paper is organized as follows. In section 2, we start by introducing our notations and provide a summary of differential geometry and kinematics for surfaces and shells. In section 3, we introduce the general governing equations: the linear momentum balance in three dimensions and resulting generic stress resultants and torque balance on a surface; the three-dimensional constitutive equation of the cell cortex, and the mass balance with turnover. Follows in section 4 the core of the theory, with the dimensional reduction of constitutive equations by asymptotic expansion, which we perform successively in the membrane and inextensional limits, based on preliminary scaling arguments. Then we introduce a variational formulation of our problem in section 5, and we derive the kinematic evolution equation for the thickness in section 6. We summarize the main equations of the model in section 7 and detail the numerical scheme employed to solve 2-dimensional thin sheet equations in section 8. We illustrate numerical results for osmotic shocks and cell division in section 9 and discuss in the light of future perspectives to extend our theoretical and numerical approaches in the Dicussion, section 10. The Appendix contains more detailed calculations and numerical values used for simulations.

Differential geometry and kinematics of surfaces and shells
Throughout the text, we use the Einstein summation convention: when an index variable appears twice in an expression, it implies summation of that variable over all values of the index. 1 Latin indices take the value 1, 2, 3, and are used in variables describing the ambient, three-dimensional, Euclidean space. Greek indices take the value 1, 2, and are used in variables in the embedded, two-dimensional, surface. Therefore, = 1 1 + 2 2 + 3 3 and = 1 1 + 2 2 . Bold symbols represent vectors. For a more complete text about differential geometry of surfaces and tensor calculus, we refer the reader to [52]; and to [53] for a concise introduction to tensor analysis. Figure 2: Parametric description of a surface patch: each set of coordinates ( 1 , 2 ) in the parametric domain  is mapped to a unique a material point on the surface  Let  be a surface embedded in ℝ 3 (e.g. the middle surface of the shell). We describe surfaces parametrically, that is, we choose a set of coordinate variables ( 1 , 2 ) in a parametric domain  ∈ ℝ 2 and let each of the ambient coordinates 1 , 2 , 3 of a material point ∈  depend on : = ( ). This corresponds to a mapping of the point ( 1 , 2 ) in the parametric domain  to a point on the material domain , see Fig. 2. This mapping is supposed to be an infinitely differentiable homeomorphism.

Differential geometry of surfaces
are curvilinear coordinates on the surface  and are non-orthogonal in general. We must, therefore, distinguish between covariant and contravariant components of tensors. At every point in the surface of the patch, the parametrization defines in-plane (covariant) tangent vectors the comma denotes the partial derivative with respect to a variable. For tangential vectors to form a covariant vector basis , we furthermore suppose them linearly independent at any point of the surface. The contravariant basis vectors are the reciprocals of the covariant basis, and are defined by the relation ⋅ = , where is the Kronecker delta. Covariant and contravariant metric tensors of the surface are defined from basis vectors as = ⋅ , and = ⋅ , respectively. These tensors allow one to generalize the notion of dot product on curved surfaces and may also be employed in lowering and raising index of tensors, for instance, = = . The determinant of the covariant metric tensor is = det( ) = 11 22 − 12 21 . The unit normal of the surface  at any point is defined by While the normal is a unit vector, the same is not true for the tangent vectors, which generally are neither of unit length nor perpendicular to each other. The second fundamental form of the surface is called the curvature tensor, defined as ≡ − ⋅ , , and is symmetric = = − ⋅ , . It quantifies the extent of local curvature of the surface in the embedding three-dimensional space.
In a curved manifold, the coordinate system is position dependent as well as the tensors living on this manifold. The covariant derivative accounts for the position dependence of the coordinate system, and extends the notion of derivatives to curved surfaces, producing tensors out of tensors. The resulting tensor is one covariant order greater than the original tensor. Covariant differentiation works differently for covariant and contravariant components and takes the following form, where we introduced the Christoffel symbols, defined by Γ ≡ ⋅ , = 1 2 ( , + , − , ) 2 . When acting on invariants 3 , the covariant derivative is equivalent to the usual partial derivative, in particular | = , . Particular covariant derivatives are the ones of the basis vectors, which are known as Gauss and Weingarten's relations, given respectively by | = , and , = − .

Differential geometry of shells
Thin shells are usually described by a middle surface (or midsurface)  and a thickness [28,27]. The middle surface is the two-dimensional surface that is equidistant, at any of its point, to the upper and lower shell boundaries. The Euclidean space in which the surface is embedded is generally called the ambient space. We suppose that the midsurface is everywhere continuous, infinitely differentiable and of arbitrary shape. In this text, we follow closely the notation of [28], to which we refer the reader for a more in-depth discussion.
The midsurface is described by a position vector ( , ) relative to an arbitrary origin through generalized coordinates . We introduce a third coordinate 4 2 Note that this definition of the Christoffel symbols suppose a torsion-free manifold (Levi-Civita connection) 3 Objects that can be constructed without reference to any coordinate systems are called invariants or geometric. Scalar fields, vectors, and tensors are examples of geometric objects 4 Derivatives with respect to the normal coordinate will be denoted as ∕ ≡ ,3 .
are the two main invariants of the curvature tensor, called the mean and Gaussian curvature, respectively. Finally, the contravariant metric tensor = is such that = .

Governing equations
The following development of the governing equations is standard and their derivation may be found in more detail in classical books on elastic shell theory [27,28]. We adopt mostly the notations employed by [28], to which we refer for more details.
At the scale of cells, the Reynolds number is very small, and we are interested in deformations at low frequencies, such that inertial and convective terms in the Navier-Stokes equation are negligible [54]. Therefore, all equations below will be written neglecting inertia. We assume that the distribution of stress throughout the sheet is characterized by the stress tensor , in equilibrium with bulk forces , described by the inertialess Cauchy's equation of motion div + = 0.
The Voss-Weyl formula [52] expresses the divergence of a contravariant tensor field without a reference to the Christoffel symbols, | = 1 √ ( √ ), when applied to the inertialess Cauchy's equation of motion it yields where is the Cauchy stress tensor and are the contravariant components of volume forces. These stresses must be expressed both per unit area of a single reference surface and relative to base vectors that do not vary across the thickness of the sheet. The natural choice for this purpose is the midsurface and its intrinsic base vectors . If we use = , 3 = , the equilibrium equations can be rewritten in the following form [28] are tangential and transverse stress tensor components referred to the midsurface. Notice that is not symmetric. Moreover, because of the symmetry of , it follows that 3 = 3 .

The cell cortex: bulk equations 3D constitutive relations
The cell cortex is often described as an active polymeric gel [55]. Polymer gels are cross-linked networks formed by linear or branched polymers. In the cell cortex, actin polymers are furthermore permanently renewed by polymerization/depolymerization, which releases elastic stresses over timescales of dozens of seconds [48,4]. The cell cortex is therefore often described as a Maxwell viscoelastic material, responding elastically on short timescales and viscously on longer times, [11]. The cell cortex differs from physical gels in which it is an active material: it is out-of-equilibrium, continuously consuming energy in the form of ATP. Here, we work in the viscous regime of cortical behavior and start from viscous active nematic gel equations [11,56,12], which consider the stress to be the sum of a passive viscous term, and an active stress exerted by myosin molecular motors in the network = + . The viscous response is modeled as a Newtonian fluid where is the hydrostatic pressure, is the shear viscosity of the actomyosin gel, and is the three-dimensional strain rate tensor, and ( , , ) is the Lagrangian fluid velocity of the shell. The cell cortex generates internal stress from the action of the motor protein myosin. Myosin motors can walk along actin filaments driven by the hydrolysis of ATP. When assembled into minifilaments and associated with several actin filaments, they induce internal stresses in the actin network (see Fig. 4). While the detailed microscopic processes for the generation and renewal of motor-generated stress are complex [57], they lead most generally to contractile stress in actomyosin networks [58,55], an effect which is akin, in first approximation, to a negative pressure. The contractile activity is controlled spatiotemporally by cells by actin myosin Figure 4: Cortex as an active gel. The myosin transforms the energy of ATP hydrolysis into mechanical work.
local phosphorylations of myosin chains through various biochemical pathways, such as the RhoA pathway [59]. As a result of motor activity or viscous flow, actin filaments may also acquire a structural organization within the cortex [60,61], that is generally measured by a director field (see Fig. 5(b)) and described with a polar or nematic 5 order parameter [11,12]. While actin filaments turn over fast, the maintenance of a local orientation may be linked, experimentally, to crosslinking [62] or filament-guided de novo assembly [63]. The active gel theory predicts an active stress contribution of the general form where Δ is the gain in chemical free energy for one ATP molecule hydrolysis,̃ > 0 is a measure of the local motor contractile activity [11]. At this point, we do not specify any particular form of the tensor , which can represent nematic or polar order 6 . Note that we neglect here the contribution of the conjugate field, a classical passive term in polar or nematic hydrodynamic theories, that reflect the restoring forces deriving from elastic energy penalizing gradients of filament orientations (the Frank-Oseen energy) [10].
The bulk constitutive relation is written in its contravariant form where we introduced, to simplify notations, ≡̃ Δ : the magnitude of the actively generated stress in the cortex, which has to be positive to correspond to contractile stress. Note that we count positively the active stress and consider a positive phenomenological parameter̃ ≥ 0, taking a different sign convention than in [11].
Local mass balance The actomyosin cortex is under permanent turnover. The polymerization of new actin monomers happens mostly at the surface in contact with the plasma membrane at = ∕2, creating a flux of actin inwards the cell at a velocity = , where is the rate of polymerization and the typical size of a monomer. The depolymerization of the actin may happen through different mechanisms [48]. 5 In contrast to polar order, the nematic phase is described by the breaking of rotational symmetry, but not translational symmetry. In the nematic phase, the molecules are randomly polarized, but their long axes are oriented on the average along a particular direction specified by a unit vector , see Fig. 5 6 In the nematic case, one can further divide the contractile stress in an isotropic and deviatoric components. This is done by simply decomposing such that The local orientation of actin filaments is measured by a traceless nematic tensor̃ = ⟨ − 1 3 ⟩ , where <> represents a local (mesoscopic) average over the orientation of the filaments in a elementary mesoscopic volume [10,64], and the orientation of each filament being characterized by a unit vector , see Fig. 5 Here we suppose that depolymerization is a bulk process that essentially happens at each free filament end, where monomers depolymerize at a rate , see Fig. 5(a). Supposing actomyosin density uniform across the membrane, bulk mass balance reads where the delta function constrains polymerization from the upper shell surface boundary = ∕2.

Balance of stress and torque resultants
The fundamental dynamic quantities of interest for a thin shell theory are stress resultants and bending moments, which are integrals of the stress and torque across the sheet. Therefore, we integrate the equations (11) with respect to through the thickness of the shell and obtain where are the stress resultant tensor and transverse shear stress are the components of the stress applied in the direction to the sheet outer surface, measured per unit area of the midsurface, and Λ ± = (1 + 1 4 ± , , ) 1∕2 . The outward unit normal vectors are ± = ± ± , with ± = −Λ −1 ± 1 2 , and ± 3 = ±Λ −1 ± . To obtain the equation for the torque resultant we multiply (11a) by , then integrate across the sheet, where d is the bending moment tensor, and  = ∫ ∕2 − ∕2 ℎ d + 2 ( + − − ) is the applied moment vector. Next, we then introduce the symmetric effective membrane stress tensor and the symmetric effective moment tensor as Then, we eliminate the shear stress resultant using equation Eq. (20), and substitute (21) and (22) into (19) to finally obtain [27] | where =̂ +  , and 3 =̂ 3 −  | . It is worth noting that the balance equations (23) are independent of the constitutive law for the shell: they represent universal force balance equations, valid for any 2-dimensional medium equipped with tension and bending moments, and simply expresses the conservation of momentum. Such decomposition is classical in the theory of elastic shells [65,27], and was also used to model the dynamics of viscous sheets [36]. Moreover, since inertia has been neglected, time does not explicitly appear in the equilibrium equations or boundary conditions.

Thin sheet constitutive equations
Because we are deriving equations of a shell theory based on stress resultants (Eq. (21)(22), constitutive equations for these stress resultants have to be formulated. These constitutive relations represent the stress field integrated over the shell thickness. Its formulation is not trivial since one must carefully consider the phenomenological three-dimensional constitutive equations for the visco-active material and for mass balance, introducing the active turnover effects into the two-dimensional constitutive relations.
We now wish to exploit the slenderness of the cortex to obtain approximated "thin sheet" equations. Such an approach is valid when the slenderness parameter = 0 ∕ ≪ 1. Then, we use asymptotic expansion in the small parameter and in to express the effective stress and bending moment in terms of the midsurface velocity ≡ ( = 0, ) = + 3 , in the local basis { , }. More specifically, the constitutive equations will depend on the two tensors that describe the rate of deformation of the midsurface: the strain rate tensor Δ , which is one half of the time derivative of the metric tensor 2Δ ≡ , and the rate of change of the curvature tensor Ω , which is defined as the time derivative of the curvature tensor Ω ≡ . These tensors can be explicitly written in terms of midsurface velocity as [36] and The deformation rate of the middle surface is described completely by the pair of symmetric tensors Δ and Ω . It is important to note that the rate of tangential stretching depends not only on the rate of surface strain ( | + | )∕2, but also on surface shape changes in the presence of curvature − 3 .

Non-dimensionalization and scaling analysis
The effective shear viscosity of the cortex is supposed to be fixed, and is a characteristic stress. Tangential coordinates scale as = ̂ , while the normal coordinate is = ̂ . Next, we define the dimensionless curvatureŝ = ,̂ = ,̂ = 2 . Then, introduce the dimensionless volume element ℎ ≡ √ ∕ = 1−2̂ ̂ +̂ ̂ 2 2 , the dimensionless shift tensor = − ̂ ̂ , and we approximate the contravariant metric tensor as, . In theories of thin sheets, one often distinguishes between two main deformation modes that are highly dependent on the load applied to the system 7 . In the membrane (or stretching-dominated) regime, the system responds mainly through stretching such that bending moments can be neglected. While in the inextensional (or bending-dominated) regime, the dissipation issued from bending deformation rates becomes significant. These two regimes are described by different scalings for velocity and stress components [34,66,32].
In the membrane limit, the velocities and stresses in a viscous sheet scale as [31] ∼ , ∼ .
While in the inextensional (bending) regime, the appropriate scaling is [33,36] In order to cover both deformation modes within our thin sheet model, we use the method of matched asymptotic expansions [49], which consists of finding approximate solutions separately in each limiting case, and then combining these approximations, giving an approximate constitutive relation that is valid for the whole range of intermediate behaviors.

Asymptotic expansion: membrane limit
We expand the variables and into a double power series in = 0 ∕ , and in the dimensionless normal coordinatê = ∕ 0 . The scaling (26) indicate the appropriate expansion in the membrane limit where ( ) and ( ) are dimensionless functions of the curvilinear coordinates 1 and 2 . We postulate that the active stress scales as the pressure, and the turnover scales as velocities (modulo a length for the depolymerization rate), which gives the following scalings for active parameters: the active stress = ̂ , the polymerization velocity = ̂ , and the depolymerization rate = ̂ . In the membrane limit, the bending moment ∼ ( 2 ), and at lowest order, = .
If the Lagrangian fluid velocity is decomposed in the local basis ( , , ) = ( , , ) + 3 ( , , ) , it follows that the strain rate tensor can be written as [28,30] 7 In the context of elastic shells, it is well established that the deformation modes also depend on the geometry and boundary conditions of the structure. For detailed description of the asymptotic behaviors we refer the reader to [30] One may remark that, at this point, no assumption was made on the kinematic, and our shell model allows for transverse shear and transverse normal deformations. Now, after applying (28) in Eq. (30), and recalling that = − ̂ ̂ , we obtain the strain rate tensor Next, we substitute the expansion (28) for the pressure, use (31) for the strain rate , and use the approximation for the contravariant metric tensor ≈ + 2 ̂ + ( 2 ). Then the leading term of the stress resultant is where we introduced the dimensionless strain rate tensor Δ (00) = 1 2 ( (00) | + (00) | ) −̂ (00) 3 , which relates to the dimensional strain rate tensor as Δ = Δ (00) .
To obtain a close relation in terms of velocities referred to the midsurface (00) , we need to find relations for the remaining coefficient (00) . These relations are found using the expansion (28) into the mass balance equation (18), and boundary condition (Eq. (102) in App. C); and then equate to zero the factors in the expansion proportional to the same power in and̂ until the leading term (00) can be fully identified. This will lead to a set of algebraic equations that can be solved sequentially.
Next, from div = , , we rewrite the mass balance (18) as where ℎ is given by (7). Then, we introduce the expansions (28), Anticipating that we only needed the lowest order terms, the expansion was truncated to ( ). The zeroth order term, ( ̂ ) 0 , is achieved when = 0, and = 0, and gives From (33), we can further obtain (00) = −2 Δ (00) − 2̂ + 2̂ By substituting Eq. (37) into Eq. (32), and rewriting the terms in its dimensional form we obtain the symmetric stress resultant where  = 1 4 ( + ) + 1 2 ,  = − 33 . We find therefore that active tension contributions are generated here not only by contractile motor activity  , as proposed earlier in axisymmetric geometry [38,7] or in general three dimensional settings [39], but also through active turnover 2 ( − ∕ ). This new active tension term due to turnover, which was overlooked so far, can lead to either a contractile or extensile stress and constitutes an important theoretical prediction of our manuscript.

Asymptotic expansion: Inextensional limit
The appropriate scaling for the inextensional regime (27) gives the following expansion and as in the membrane regime, we hypothesize that active sources of stress scale similarly to the sheet dynamic variables: = 2̂ , = 3̂ , and = 3̂ . We then use a similar approach to obtain leadingorder expressions for the dimensional (symmetric) stress and bending resultants , . For the sake of the presentation, we defer detailed calculations to Appendix A.
We find for the tension resultant where Note that in this regime the active turnover stress appear only in the terms proportional to 3 in the membrane stress, and in the bending moment.

Composite equations
The composite equations for and are obtained by gathering together the different terms appearing in the membrane and inextensional limits. The approximate constitutive relation obtained by "matching" the terms is valid for the whole range of intermediate behaviors. The ensuing complete system is then, In the absence of active terms, our constitutive relations are identical to the ones obtained in [36] for a viscous sheet without couplings between bending and stretching. Moreover, we note that the absence of transverse shear deformation at the lowest order in the asymptotic expansion induces a Kirchhoff-Love type of deformation 8 . For an initially curved elastic shell, Kirchhoff-Love deformations lead to a theory generally named after Koiter [67]. This can be readily seen since the passive parts of the resultants (42) and (43) are, modulo a time derivative, similar to the incompressible limit of Koiter's constitutive equations for elastic shells [67], as one may expect from the Stokes-Rayleigh analogy 9 .

Variational formulation: weak form and Rayleigh potential
The weak form is obtained by multiplying the balance equation (23) by a test function, then integrating over the region in space occupied by the shell, where ( 1 , 2 ), and 3 ( 1 , 2 ) are a set of "virtual velocities" (test functions). Using integration by parts, the divergence theorem for curved surfaces, and neglecting the work terms from the surface edges, the variational (weak) form is written as [27] where , and are defined in (21), and (22), and where Δ = 1 2 ( | + | ) − 3 , and Ω = 3, | + | + | + | − 3 . The left hand side of the weak formulation (45) might be interpreted as internal virtual power, since we identified Δ and Ω as measures of strain rate increments. The terms in the right hand side of (45) can be interpreted as the virtual power of external loads.
The Rayleigh potential Starting from the weak form (45), if the set of virtual velocities are chosen to be the real velocities = , one can identify the resulting variational function as a Rayleigh potential [68,69], , which is the sum of the (passive) power dissipated by viscous forces , and the power input of active terms .
In the absence of external forces, the Rayleigh potential is The Rayleigh potential plays a similar role for viscous problems as the energy potential in elasticity, and allows to calculate the midsurface velocity as its minimizer The midsurface geometry may then be deduced from the current velocity as = . In a sense, this simple variational form expresses here simply the fact that viscous flows minimize dissipation [54]. A generalization of such approach for dissipative systems, to include other sources of dissipation but also reactive forces deriving from a free energy, is often named Onsager variational principle in soft-matter physics [70,71].

Evolution of the sheet thickness
To obtain the equation of evolution for the sheet thickness, one may integrate the equation (34), as performed in [36]. Alternatively, we derive here mass balance directly at the level of the shell. Consider a small volume of the sheet d = ∫ In the left hand side, we use Leibniz's integral rule to bring the derivative inside the integral, and since the integration limits ± ∕2 are time-dependent, this will give terms in ∕ . Then, using the fact that = (det ) = det , it follows from (24) Given 2 ∼ ( 2 ) ≪ 1, the equation for the time evolution of the thickness becomes, at first order in For an incompressible passive fluid, the equivalent equation would be simply = − Δ , which states that rate of change of the thickness is proportional to minus the shell stretching rate Δ , as found in [33,36]

Summary of the main equations
Solve for the midsurface velocity that minimizes the Rayleighian, with where the thin shell constitutive equations are The midsurface velocity is implicitly written in terms of rate of deformation measures, namely the strain rate tensor and rate of change of the curvature tensor The geometry is then solved in a Lagrangian manner, the midsurface geometry position is deduced from the current velocity and the thickness spatio-temporal evolution is described by the following equation

Preliminaries
We discuss here some aspects regarding the finite-element solution scheme used for the numerical simulations. Note, however, that computational aspects of the viscous thin shell theory are not the primary purpose of the present work. The proposed implementation will therefore be limited to simple geometrical settings which excludes non-manifold surfaces for instance. For more details on the practical numerical implementation aspects, we refer to the documentation of the source code released with the paper [50].
For viscous membranes, a few numerical implementations have been recently proposed based either using finite element models [72,39] or discrete geometric approaches (computer graphics) [73,25]. Yet, no numerical work has been dedicated to implementing viscous thin shells in generalized coordinates to our best knowledge. Let us first point out that thin shell theories based on the Kirchhoff-Love kinematics are notoriously difficult to implement due to the system being of fourth-order type, which requires 1 -continuity of the approximating finite element discretization fields to avoid "kinks," leading to infinite bending energies. Since this cannot be achieved with standard  0 finite elements, one has to resort to either dedicated or complex finite elements [74,75,30] or isogeometric analysis approaches [76]. To circumvent the difficulty of handling  1 -conforming elements, one may also relax the underlying thin shell hypothesis and solve an equivalent shearable shell model [77] which requires only  0 continuity. The use of such a shearable (thick shell) instead of unshearable (thin shell) kinematics must, however, consistently converge to the latter in the thin shell limit → 0. In such thick (or Naghdi-type) shell models, the shell cross-section does not remain normal to the shell surface any more. Additional degrees of freedom are introduced in the form of inextensible directors. These inextensible directors are allowed to deviate from the surface normal during the deformation, hence accounting for transversal shearing of the shell across the thickness. More simple  0 elements can then be considered in the case of thick shell kinematics for the both displacement/velocity and director fields. However, shear locking issues are prone to arise in the thin shell limit, requiring the use of dedicated numerical techniques such as mixed methods to circumvent the locking issue [78,79].
Our simulations were conducted along these lines, by relaxing the thin shell hypothesis in the numerical implementation by considering a director field which, contrary to , is not necessarily normal to the shell midsurface, see Fig. 6. Strain measure expressions need then be adapted with instead of . Since the shell cortex is, in practice, extremely thin, we expect to recover the theoretical thin shell limit derived above by adding a penalization for the deviation between the normal and the shell director in the dissipation potential (see Eq. (63)).

Discrete shell geometry
The discrete shell reference surface, denoted by , is assumed to be made of an assembly of flat triangular elements and 0 = 0 denotes the reference shell director which we assume to initially coincide with the reference shell normal. The reference metric tensor and curvature tensors are defined respectively as 0 = Reference configuration Current configuration Following [79], we assume the shell director to be unstretchable (‖ ‖ = 1) so that we parametrize the unit director by two angles ( 1 , 2 ) (see Fig. 6 right), giving rise a so-called five-parameter shell model: [80] ( 1 , 2 ) = (cos 1 sin 2 , − sin 1 , cos 1 cos 2 ). (59) Note that this simple parametrization limits our implementation to non-closed manifold surfaces. We implement an updated Lagrangian formulation, in which the discrete equations are formulated in the current configuration, which is assumed to be the new reference configuration. In this setting, we replace the definitions of the strain rate tensor Δ = 1 2 and rate of change of curvature tensor Ω = with their linearized expressions for the rate of deformation [79]: and where ≡̇ = is the rate of change of director, and = . In addition, we must also consider the rate of shear strain, When = 0, we recover that = ( 3, + ) , then substituting into Eq. (61) we should retrieve exactly Eq. (25), which corresponds to the unshearable shell dynamics 10 .
As a consequence, the new Rayleigh potential density is modified to add a novel dissipation contribution penalizing shear: where ,pen is a shear penalization constant, that we choose sufficiently large to ensure ∼ everywhere on the midsurface.

Finite-element discretization
At each time step, we must solve for the shell midsurface velocity and the new director . In fact, we consider the rate of changė of the -parameters to be the main unknowns from which we also have = ⎡ ⎢ ⎢ ⎣ − sin 01 sin 02 cos 01 cos 02 − cos 01 0 − sin 01 cos 02 − cos 01 sin 02 As a result, we have a system with 5 degrees of freedom to solve at each time step: three velocity components ( 1 , 2 , 3 ) = and two angular rates of change (̇ 1 ,̇ 2 ) =̇ . In the present work, we do not follow the MITC discretization discussed in [79] to remove shear locking but rather use a much simpler finite-element discretization, namely: • a quadratic continuous ( 2 -Lagrange) interpolation for the velocity field • a linear interpolation with continuity enforced at the triangle mid-sides (Crouzeix-Raviart element) for the angular rateṡ Such a discretization choice is extremely simple and seems free of any shear locking effects as shown in [81] for nonlinear shells and [82] in the case of plate and shell limit analysis theory. 10 We use the fact that , = ( | − 3 ) + ( 3, + ) , which comes from the decomposition of the velocity in the local frame and the use of the Gauss-Weingarten relations

System evolution
When considering the system evolution between time and +1 , the resulting linear system for the current velocity and rate of change of director field is solved when considering the shell thickness to be fixed to its previous value i.e. = +1 ,̇ +1 = arg min ,̇ ( ,̇ ; ) (65) where is the total Rayleighian of the system, of surface density given by (63). Here, we assume no external loading on the shell, such that all shell deformations are powered by internally generated active stresses. Note that, due to the linear expressions of the strain rates (60), (61), (62) and the quadratic form of the dissipation potential , the resulting system to solve is, in fact, linear [83,7]. We then update the geometry in an asynchronous and Lagrangian manner, by solving the discretized thickness evolution equation and updating the midsurface position and director using (59) as follows: The full finite-element implementation is done within the FEniCS software package [84,51] and the meshes are produced with the Gmsh generator [85]. In the geometry updating process, we also regularly call some remeshing procedure using the MMG platform [86]. The Python code used for implementation in FEniCS, as well as the two following numerical illustrations may be found here [50]. Considering the symmetry of the numerical examples presented in Sec. 9, the finite element model of the cell is modelled as one-eighth of a sphere.

Volume constraint
In some numerical illustrations, such as cell division, Sec. 9.2, we want to ensure that the total volume of the cell is conserved during the division. To do so, we enforce the following constraint, Enforcement of this constraint is achieved through the introduction of a single scalar Lagrange-multiplier which can be interpreted as the cell hydrostatic pressure .
Note that the previous approach can also be easily extended to the case where the volume change rate is imposed to a given constanṫ imp (71) by considering this new Lagrangian instead:

Numerical illustration: osmotic shocks and cell division
This section presents some biologically relevant model problems illustrating the partitioning of the sheet deformation between stretching and bending dominated modes. In particular, cell division and hypo-osmotic shocks are mainly dominated by tensions, whereas fast hyper-osmotic shocks will be bending-dominated.
In the following, for the sake of simplicity we neglect local variations in the filament orientation within the cortex. We furthermore suppose that filaments remain parallel to the cell surface, which is a good approximation for the true organisation of actin in the cell cortex, as observed experimentally [87,88]. Since we neglect tangential gradients of orientation, filaments are supposed isotropically oriented along the midsurface (see Fig. 5). To account for the isotropic organization of filaments along the midplane, we assume a purely nematic form of the active stress tensor in 3D, with constant diagonal coefficients in the local midsurface basis ( , ) [89,12,7 In all generality may also have an isotropic contribution. However such isotropic term would appear only in active bending moments, but not in membrane resultants. We will see that in our two examples below (osmotic shock and cell division) that active bending moments play a negligible role in the shape dynamics, and we therefore ignored such contribution. To explicitly consider filament reorganization one may have to write down and expand a dynamic equation for the spatiotemporal evolution of this nematic tensor and its coupling to velocity gradients [11,12,9,90,91], which is out of the scope of this manuscript.

Osmotic shock
Here, we mimic the response of a suspended cell to rapid osmotic shocks. Considering a cell in osmotic equilibrium with an external medium, a sudden change in solute concentration creates an osmotic imbalance that will generate a rapid flux of water across the cell membrane, changing its volume until a novel equilibration of concentrations is reached. If the solute concentration in the medium is increased, the cell is subjected to a hyperosmotic shock, causing an outward flux of water and cell shrinkage. In contrast, a decrease in solute concentration leads to an inward flux in the cell called hypo-osmotic shock, leading to its expansion. For the sake of simplicity, we don't intend here to solve for a realistic pump-leak model of cell volume control, involving electroosmotic effects [92,93,94], but we impose instead directly the variation of a target volume ( ) using the constraint (72) (see Figs. 7a and 8).
We first consider an hypo-osmotic shock, where a solution may be calculated analytically. We suppose that the cell's response is mainly in the membrane regime, and that the volume can be externally controlled and is given by ( ). Because the deformation is spherically symmetric, the time evolution for the radius of the expanding cell reduces to ( ) = 3 √ 3 ( ) 4 . In spherical coordinates, the metric and curvature tensors are = diag( 2 , 2 sin 2 ), and = diag(− , − sin 2 ), respectively. Moreover, since Δ = 1 2 , we have for isotropic deformations Δ = diag( ̇ , ̇ sin 2 ). The thickness evolution, Eq. (50), can be readily written aṡ The first-order differential equation (74), subjected to the initial condition ( = 0) = 0 can be analytically solved for a given ( ), consequently, given ( ). The steady-state thickness at a given cell radius is given by 0 = − . It slightly differs from the value ∕ value defined in [7,12], which is valid only in the limit of a flat plate → ∞. It is corrected here by a curvature term − in the thickness evolution, that comes from the asymmetry of the polymerization and the curved geometry of the sheet.
Such simplified system can be used to validate the FEM approach: we impose a sigmoid volume variation ( ) and compare the analytical solution of (74) for the thickness ( ) with the solution obtained through finite elements under the same set of parameters, see Fig. 7.

Young-Laplace equation for thin shells
In a spherical symmetric system, assuming that tangential tension gradients and bending moments may be neglected, an analog of the classic Young-Laplace equation can be retrieved for thin shells from the out-of-plane equation (23b), where is the internal pressure that maintains the cell volume constant (external pressure is supposed zero), and the tension supposed isotropic and homogeneous along the membrane. It is interesting to note the geometric correction to Laplace's law in ∕ = ( ) in front of the pressure, that is usually neglected for very thin shells or membranes.
In the approximation of slow volume change (̇ = 0,̇ = 0) the above coupled equations at steady-state reduce to where f lat 0 ≡ is the stationary thickness in the limit of a flat shell → ∞ with no flow. One can remark that the typical Laplace's formula used in spherical geometry = 2 f lat 0 2 [12,7] is corrected by a term of order ∕ = ( ) proportional to the turnover rate .
Hyperosmotic shock To probe our model in situations where bending dissipation may become non-negligible, we study the effects of a rapid hyperosmotic shock on a cell. Here, turnover acts as a stabilizing agent: if the volume variation is slower than the typical turnover timescale, the cell cortex will have the time to rearrange and adapt the cell's shape to its new volume without buckling. We thus apply a sudden drop in volume that we model by a sigmoid function, as illustrated in Fig. 8. As the cortex does not have time to depolymerize homogeneously, folding induces spatial variation in the thickness correlated to local variations in curvature. In our simulations shown in Fig. 8, we find numerically that active bending terms play a negligible role, contributing to less than 0.1% of the total power input. However, passive bending terms are not negligible. To compare the relative contributions to the dissipation of bending and stretching rates, we define the following non-dimensional parameter which measures the relative weight of bending dissipation  = ∬   1 3 3 Ω Ω d to the total viscous dissipation in the shell  =  +  , where  = ∬  4  Δ Δ d is the dissipation due to stretching. The time-course of the parameter  is plotted in Fig. 9. Two distinctive peaks indicate a significant contribution of bending to dissipation, which may reach up to 60% of the total. Interestingly, these two peaks, indicated by (I) and (III), do not correspond to the most folded cell shape, which rather lies at the minimum of the curve between these two maxima, as indicated by (II). This is because, for a thin viscous sheet, dissipation increases with the rate of curvature change, but not with the curvature itself, like one may be used to in elasticity. Therefore, maintaining highly curved cell shapes does not cost energy, but dissipation shall oppose rapid shape changes.

Cell division
Cell division encompasses two main processes in biology: mitosis, which corresponds to the even splitting of genetic material between the two future daughter cells, and cytokinesis, whereby the mother cell cytoplasm is physically cleaved in two [95]. In the following, we focus on the cell's shape deformation during Figure 9: Plot of the relative contribution of the bending dissipation to the total dissipation. The two peaks (I and III) represent the instants in time for which the rate of curvature change is the highest. Even though at intermediary time II the cell shape is overall much more curved, bending dissipation is the lowest.
cytokinesis. In animal cells, cytokinesis is achieved through the ingression of a circumferential actomyosin furrow located at the equator of the cell, while cortical tension at cell poles is resisting the deformation imposed by increasing cell pressure [96,7]. Various mechanical models of cell deformation over cytokinesis have been proposed in the last decades, with different assumptions about the mechanical nature of the cortical material: contractile [97], elastic [98], viscous [99] or viscoelastic [100]. Other approaches have considered the lipid bilayer mechanics only and attributed most of the mechanical dissipation to cytoplasmic and surrounding medium fluids [101,102]. In some of these models, a contractile ring structure of different nature is supposed to provide the constricting force at the cell equator [98,100], while in others, the furrowing is the natural result of higher but continuous activity in the equatorial region, that may trigger possible structural rearrangements of the continuous cortical surface [97,99]. In agreement with early hypotheses [103], recent active-gel models of cytokinesis [7] and experiments [96,61] all point towards the continuum character of the cortex, that may flow in a Marangoni-like process toward the equator as a result of contractile gradients and can reorganize structurally along the cell's surface in response to extensile and compressive forces. In fact, the hypothesis of surface tension gradients driving division is ancient and goes back to early observations of dividing marine eggs [104,105,106] and to the similarity of their shape to one of the oliveoil droplets in response to surfactants [97]. A gradient of contractility centered at the equator has indeed been measured from the cell poles to the cell equator in the form of a Gaussian-like spatial distribution of active RhoA, a protein upstream of actomyosin contractility [107].
Following [7], we start from a spherical cell of radius 0 and mimic the RhoA gradient by imposing a Gaussian distribution of active stress along a fixed ambient coordinate direction centered at the equator, where 0 is the basal level of stress, that creates an active tension at poles, and is called equatorial overactivity that expands over a typical width . In the initial state, the cell is a sphere of 0 and uniform thickness = ∕ , which is subjected to an uniform contractile activity 0 . A typical deformation timescale in this system is the ratio of the viscosity to the contractile stress ≡ ∕ . Typical numerical values for all these parameters are extracted from experimental measurements referenced in the literature and may be found in Table 1. Following the rise of a Gaussian over-activity, the equatorial region of the cell becomes more contractile and generates cortical flows toward the equator. This triggers a local increase in the thickness, which is partially counterbalanced by depolymerization, as found in the mass balance Eq. (34). Material turnover is indeed the essential factor ensuring the overall stability of the thickness evolution in presence of contractility or thickness gradients. The accumulation of matter in the furrow, here represented by the increase in thickness, increases further the contractile stress around the equator 11 . The resulting contractile ring then pinches the cell, leading to its constriction. The time course for such dynamics is shown on Fig. 10. The spatial distribution of contractility is assumed here to be time-independent, but it may slightly vary with cell shape changes [107]. The spatiotemporal evolution of the thickness and of the midsurface velocity field are illustrated at corresponding timepoints in Fig. 11.  As for the osmotic shocks, we find numerically that active bending plays a fully negligible role. But, in contrast to the hypoosmotic shock, bending dissipation remains here also largely negligible during the timecourse of cell division (not shown), and the cell adopts very curved shapes in tje final stages of cytokinesis with almost no associated dissipation. Cytokinesis may therefore be well described by a simple membrane theory, where tension is the dominating phenomena, which justifies a posteriori the approach of [7], where an axisymmetric active viscous membrane model was employed to simulate cytokinesis. As originally predicted in [7], we show that full constriction of the furrow, leading to the completion of cell division, depends on the level of equatorial over-activity ∕ 0 , presenting a threshold , below which the full constriction fails. We plot in Fig. 12b, the bifurcation diagram for physiological values of the parameters. Each point in the figure represents the final radius following the dynamics as shown in Fig. 12a, with variable signal over-activity.
Even though the obtained results in [7] were qualitatively correct, capturing well the shape changes for axisymmetric cell division, the model neglected an essential stress contribution, generated by active material (de)-polymerization. To compare the relative contribution of the turnover to the total input power, we define 11 Simple stability analysis of a flat shell under uniform active stress shows that the thickness is stable if ≤ 4 . If the additional tension generated by turnover 2 ( − ∕ ) is neglected, the threshold for stability doubles to ≤ 8 , which shows that the turnover stress has a destabilizing effect. the following dimensionless parameter where  = ∬  − Δ d , and  = ∬  1 2  Δ d are active power contributions due to turnover and contractility, respectively. Note that we did not include the active bending term in  , as its contribution to the active power input is negligible. The stress generated by the polymerization comes from the local difference in thickness 1 − 0 created by the cortical flows with a reference thickness 0 = for a flat plate with vanishing flow. Since the variation in thickness is highly concentrated around the equator of the cell, most of the active power input due to the turnover stress will be concentrated in this region and accounts to up to 25% of total (integrated) active power input, as shown in Fig. 13. Time t/t a Local turnover stress Figure 13: Relative contribution of the (integrated) power input of turnover with respect to the total power input due to both turnover and contractility as function of time. The insets show that most of the stress generated by the polymerization is localized at the furrow.

Discussion
The lack of generic and flexible theoretical and numerical frameworks to realistically simulate the 3D shape changes of single or multiple interacting cells has been limiting our ability to characterize complex morphogenetic events, such as asymmetric cell divisions [108,109,110,111,112] and rapid multicellular reorganizations [113,114,115], which are fundamental to early embryogenesis and tissue morphogenesis. To fill this gap, we derived in this paper a viscous active thin shell theory for the cell cortex in general curvilinear coordinates, and we proposed a numerical implementation with the finite-element library FEniCS [84,51], which we hope will facilitate the spread, reuse, modification, and refinement of our model.
We have used systematic asymptotic expansions to derive leading-order constitutive equations governing the dynamics of the actomyosin cortex. The main geometrical assumption is the slenderness of the cortex, meaning that the characteristic thickness 0 must be small relative to a characteristic size of the cell, typically its curvature radius. Our main physical assumption is to describe the actomyosin cortex as a continuous medium, relying on hydrodynamics active-gel models [11,9]. Here, we chose to model the actomyosin material as a single incompressible active fluid phase following previous approaches [12,7]. However, the full rheological and dynamic complexity of actomyosin networks may not be fully encompassed by such a simple constitutive equation. The asymptotic expansion approach may provide a generic and rigorous method to derive alternative thin shell models based on refined descriptions of actomyosin rheology. A first extension would be to account explicitly for the dynamics of the nematic order of filaments. Active nematic gel theories [9,12] provide a convenient framework to describe the bidirectional coupling between nematic order and surface flows. Our derivation was done with a generic active tensor . It was then supposed nematic in the two examples to account for the approximate alignment of filaments parallel to the midsurface, while we neglected dynamic filament reorganization in this tangent plane. To obtain a proper nematic active thin viscous shell theory, one should write explicitly a bulk dynamic equation for and reduce it into an effective 2-dimensional equation defined on the midplane. Our theory may also be readily extended to account for possible transverse gradients in myosin [116], that could generate additional active torques. We assumed that actomyosin is incompressible, but this material may be better described as a permeating gel by introducing a local actin density [117,118]. Describing the cortex as a gel seems important to better assess the contribution of turnover to stress resultants, which relies here importantly on the incompressibility condition. By essence, our model was designed as purely viscous and is well-suited to describe shapes changes on timescales much larger than a typical actomyosin turnover time −1 . But it could also be extended to describe the viscoelastic response of the cortex on shorter timescales [4], although this may entail the more difficult handling of power-law rheological models [119]. All these possible extensions rely generally on bulk descriptions of actomyosin mechanics that do not derive directly from microscopic descriptions of the network. Ample numerical and theoretical effort will also be required in the field to bridge these scales quantitatively.
Our thin shell theory accounts for viscous stretching and bending resistance, making it relevant to describe rapid curvature changes. Although viscous in nature, the Stokes-Rayleigh analogy allows for interpreting our model in terms of the classic theory of elastic shells, and the equations we obtained may be classified into Koiter-like shell models. In contrast to alternative direct active surface approaches [120], the derivation of a proper thin shell theory starting from three-dimensional constitutive equations, shows its relevance for the cortex through two important predictions: • First, the balance between polymerization and depolymerization is predicted to create contractile or extensile tensions and may also oppose or promote the bending of the layer. This effect depends inherently on the shell character of the cortex and relies primarily on our hypothesis of an incompressible phase and on the fact that the cortex grows from the plasma membrane, creating an inward addition of material that may translate into longitudinal viscous flows, and hence into turnover-generated viscous stresses. It follows that active cell deformations depend not only on the contractile activity of myosin motors but also on the speed of actin (de)polymerization, which is particularly important to relate cell shape changes to its molecular regulation. Such contribution has been neglected in previous theories of active membranes [12,7,39]. Yet, using experimentally measured parameters, we find that up to 25% of the active power input may come from turnover during cytokinesis and that it largely contributes to the equatorial contractile force that divides the cell.
• Second, a proper thin shell model with bending resistance is necessary to account for rapid changes in cell curvatures, as one may observe in hyperosmotic shocks. Here, the asymptotic approach allowed us to calculate prefactors for active bending terms (43), which would have been otherwise difficult to write down directly. Yet, in the two numerical examples we chose, we found that active bending was largely negligible.
From a numerical perspective, we provided an easily modifiable Python implementation of the model within the general-purpose finite-element library FEniCS. We proposed a numerical scheme based on a penalization method that circumvents the need for 1 elements to use more simple 0 elements. In the cortex, gradients in contractile stress can lead to strong tangential motions, and convergent flows, which create excessive distortions of the computational mesh linked to the material [7] (see Fig. 10). In the present work, such significant distortions were circumvented through repeated remeshing of the computational mesh and interpolation of the fields. This solution may prove costly for larger meshes, and a more appropriate alternative using an Arbitrary Lagrangian-Eulerian (ALE) representation [121] may be used, in which one is free to decouple the movement of the Lagrangian particles from the movement of the mesh. Such an approach has been previously employed for the simulation of fluid membranes [39,122] and our description can be very easily translated into a similar ALE formulation.
Our finite-element numerical scheme may be adapted in the future to tackle more complex morphogenetic processes. In particular, one can readily envision the study of the cell division when the cell is confined to a rigid shell, as happens for instance in sea-urchin or worm eggs, such as C elegans, where such confinement leads to two tightly apposed membranes in the furrowing region [123]. During division and migration, cortical flows have been long predicted to reorient actin filaments [103,60,61], which in turn may influence cortical contractility [97,99]. Interestingly, Rayleigh variational formulations for 2D nematic surfaces are already available [90] but solving numerically for the coupled viscous-nematic dynamics remains a challenging task [91], that may require specific finite-element approximations to implement tensorial dynamics on moving boundaries [124,125]. Another possible extension of our numerical approach is the coupling of cortical mechanics with regulatory proteins lying at the plasma membrane, which may diffuse, react and be advected by cortical flows to generate complex morphogenetic patterns [126,127,128,129,38,19]. Finally, the generalization our model to multicellular systems, by solving, for instance, similar equations on non-manifold triangular meshes [8], would allow us to envision very realistic simulations of cell mechanics in early embryo development. = − Δ (20) . Then, using Eq. (92), we obtain (10) = −2 Δ (20) .

B. Numerical values
The physical parameters used in the main text are summarized in Table 1.