Stability of asymmetric cell division under confinement: A deformable cell model of cytokinesis applied to C. elegans development

Cell division during early embryogenesis has been linked to key morphogenic events such as embryo symmetry breaking and tissue patterning. It is thought that boundary conditions together with cell intrinsic cues act as a mechanical “mold”, guiding cell division to ensure these events are more robust. We present a novel computational mechanical model of cytokinesis, the final phase of cell division, to investigate how cell division is affected by mechanical and geometrical boundary conditions. The model reproduces experimentally observed furrow dynamics and predicts the volume ratio of daughter cells in asymmetric cell divisions based on the position and orientation of the mitotic spindle. We show that the orientation of confinement relative to the division axis modulates the volume ratio in asymmetric cell division and quantified the mechanical contribution of cortex mechanics, relative to the mechanical properties of the furrow ring. We apply this model to early C. elegans development, which proceeds within the confines of an eggshell, and simulate the formation of the three body axes via sequential (a)symmetric divisions up until the six cell stage. We demonstrate that spindle position and orientation alone can be used to predict the volume ratio of daughter cells during the cleavage phase of development. However, for compression perturbed egg geometries, the model predicts that the change in confinement alone is insufficient to explain experimentally observed differences in cell volume, inferring an unmodeled underlying spindle positioning mechanism. Finally, the model predicts that confinement stabilizes asymmetric cell divisions against bubble-instabilities, which can arise due to elevated mitotic cortical tension. Author summary A crucial morphogenic step during early embryonic development is symmetry breaking in the embryo. For C. elegans the formation of the three body axes can be traced back to the six cell stage, where tissue-topology is the result of symmetric and asymmetric divisions. How cell mechanical boundary conditions and cell intrinsic cues influence this process of symmetry breaking is still an open question, as currently, a quantitative mechanical description of cytokinesis in complex architectures is lacking. We developed a simple mechanical model of cell division, incorporated in an existing mechanical cortex model, to simulate cytokinesis in geometrically confined environments. Our approach was able to both capture furrow ring dynamics and predict the volume ratio of daughter cells accurately. By simulating early C. elegans development with different geometrical boundary conditions, we were able to trace back the origin of volume discrepancies between the experimental setups to a quantifiable shift in spindle positioning during cytokinesis. Finally, we showed how embryo confinement partially stabilizes bubble-instabilities that arise during asymmetric cell division during the early cleavage phase.

Early embryogenesis is characterized by a well-orchestrated sequence of cell division and 2 differentiation events. In bilateria, an important step during this process is the 3 establishment of the three main body axes: Anterior-Posterior, Dorso-Ventral and 4 Left-Right [1]. In model organisms with an invariant developmental cell fate, such as 5 the nematode Caenorhabditis elegans, the origin of these symmetry breaking steps has 6 been traced back to asymmetric cell divisions, which can be biochemical and/or physical 7 in nature [2,3]. Although the importance of physical asymmetry, such as an asymmetric 8 volume ratio of daughter cells, remains poorly understood, recent observations for P0 9 divisions in C. elegans have linked daughter cell size asymmetry to changes in temporal 10 cell cycle coordination, tissue topology, cell fate, and even overall robustness of early 11 embryogenesis [3]. 12 To control the volume of daughter cells, cell division relies on the precise positioning 13 of the spindle apparatus, as the future cleavage site aligns with the spindle 14 midplane [4][5][6]. Accumulated discrepancies during this process can lead to polyploidic 15 cells, which have major detrimental effects on tissue homeostasis and development [7][8][9]. 16 It is thus imperative for a dividing cell to have strict spatiotemporal control of its 17 spindle, as any shift in spindle position influences the volume ratio daughter cells and 18 viability [10]. Multiple spindle positioning mechanisms have been identified and 19 categorized based on the nature of the underlying cues e.g., internal vs external and 20 biochemical vs mechanical [11]. Spindle positioning in the P0 division of C. elegans for 21 example, which is a highly asymmetric both in volume and distribution of molecular 22 species, is governed by a biochemical gradient which induces a bias in the magnitude of 23 cortex-spindle force interactions towards the posterior pole [12]. For wild type C. 24 elegans, this bias produces a net spindle division plane offset, from the cellular 25 center-of-mass in the direction of the posterior pole. This offset eventually results in an 26 asymmetric division where the anterior cell is significantly larger than the posterior cell, 27 at a ±60/40 volume ratio. Disrupting biochemical cues, e.g. the distribution of PAR 28 proteins, or interfering with mechanical cortex-spindle interactions can lead to 29 symmetric P0 division [3,13]. Interestingly, C. elegans embryogenesis seems robust to 30 compression induced eggshell deformation [14]. However, in these geometrically 31 perturbed setups, changes in cellular movement [14] and enhanced volume asymmetry in 32 divisions that contribute to the LR-axis [15] have been observed in the early embryo. 33 Mechanically, the final phase of cell division, cytokinesis, is driven by a tensile 34 furrow ring consisting of cross-linked actin filaments [16]. The dynamic cross-linking of 35 tread-milling actin subunits in the ring is thought to drive tension build-up, resulting in 36 furrow ingression [17]. However, the specific contribution of actin-motor cross-linkers 37 such as non-muscle myosin-II is an ongoing debate [16]. As of yet, while key biochemical 38 mechanisms that regulate cytokinesis have been identified [2,[18][19][20][21], a full 39 spatiotemporal mechanical description of cytokinesis at cell scale is still lacking. The 40 past decade has however seen significant progress in the development of mathematical 41 and computational models that describe the physics of cytokinesis. However, most fail 42 to account for the geometrical and mechanical constraints imposed by the surrounding 43 tissue architecture and the embryonic boundary. Furthermore, few mechanical models 44 March 15, 2022 2/14 exist that are able to represent the physics of asymmetric cell division [22,23]. Yet, 45 especially in early embryonic development, it is thought that the coupling between the 46 asymmetry of cell divisions and the geometry and connectivity of the surrounding tissue 47 stabilizes robust tissue patterning [24][25][26][27][28]. At long timescales, the model of a cell as a 48 fluid droplet with surface tension and adhesive bonds provides a minimal description of 49 the mechanical response of a cell. It can be tuned with a small set of experimentally 50 accessible parameters, such as the interfacial tension and the effective cortical 51 tension [29,30]. However, droplets under tension are inherently susceptible to 52 "bubble-instabilities", where regions with locally elevated curvature, such as the small 53 daughter half in an asymmetric division, tend to collapse in favor of low curvature 54 regions due to differences in local pressure. Interestingly, such instabilities have been 55 observed in vitro in HeLa cells and fibroblasts [31], where they produce shape 56 oscillations during cytokinesis. 57 We present a novel computational model of cytokinesis based on an existing 58 Deformable Cell Model (DCM) [32][33][34][35]. In the DCM, the cell cortex is approximated as 59 a viscous shell under surface tension and cell-cell interactions are takes into account 60 through adhesion forces. A complete description of the DCM can be found in S1. For 61 this work, we made abstraction of the underlying molecular mechanisms of furrow ring 62 mechanics and used a coarse-grained representation where the furrow ring was 63 approximated by a closed loop consisting of elastic subunits, which shrink as cytokinesis 64 progresses [36], see Models section. Furrow ingression then follows as a result of elastic 65 tension build-up in the contractile ring and local cortex deformation. Ring contraction 66 progresses until cell abscission is initiated at the prescribed abscission radius. Following 67 cell abscission, membrane and cortex wound healing wrap up cell division, splitting the 68 mother cell into daughter cells [37]. C. elegans was used as our biological model system 69 due to its invariant developmental path and transparent eggshell, allowing for reliable

75
We first validate our cell division model by comparing simulation results to 76 experimental data and theoretical predictions found in literature [6,17,31,36]. We start 77 by investigating furrow ring dynamics, before focusing on spindle positioning and the 78 volume of daughter cells. Triangles colored in red are marked as furrow ring. In these region, tension is generated to drive furrow ingression. After abscission, the cell is split into two daughter halves and a new computational mesh is generated. Schematic representation of a dividing cell. Triangulated surface meshes are used to represent the viscous cortical shell with thickness t, passive mechanical properties (viscosity η) and active mechanics (contractility γ and bulk modulus K ) [33,34]. Distinction is made between interphase γ I and mitotic surface tension γ M , fixing γ I = γ M /3. The furrow ring, represented in red, virtually splits the mother by generating contractile forces per sub-unit, resulting in a net inward facing force. B: Contractile ring perimeter in function of time during cytokinesis. Larger cells spend proportionally more time in the constant constriction regime. C: Constriction rate in function of ring perimeter. A constant and linear constriction regime can be identified, in line with experimental observations. D: Total cytokinesis duration in function of mother cell radius, together with Eq. (1), normalized with the theoretical division duration. which put forward an estimate for the total duration of cytokinesis as with τ d as the base ring closure time and τ t the expected furrow ring closure time.

95
Model input parameter τ d = 263 s was estimated based on the work of Carvalho et 96 al. [36], see Models section.

97
As our simulations were able to reproduce experimentally observed scaling of τ t , see 98  For spherical cells, the relative offset of the division plane, α n = α/R 0 , with respect 113 to the center-of-mass and daughter volume is given by where V n represents volume fraction of the largest daughter cell with respect to the 115 mother cell. For regular geometrical shapes, similar relations can be derived. We properties drive the bubble-instability of an asymmetric division, i.e., any asymmetry 132 will be further increased, potentially leading to complete collapse of the smaller 133 daughter half. We postulate that a biological cell escapes this instability through a 134 non-zero bulk rigidity of the cytosolic content, which we implemented in our model by 135 means of a daughter-half bulk modulus K div . In accordance, we define the relative 136 mitotic tension compared to K div as γ n := 2γ M /(R 0 K div ), a measure that expresses We apply our mechanical model for cytokinesis to asymmetric cell divisions that take 148 place during early C. elegans development. We focus on the cleavage phase up to the six 149 cell stage, as development up to this stage is thought to be predominantly the result of 150 passive cell mechanical and cell-cell interactions, and less of active mechanical processes 151 such as cell migration [40]. We investigated both compressed and uncompressed embryo 152 shapes. Embryo compression occurs naturally in-utero, while in-vitro it generally occurs 153 as a result of sample preparation for imaging, where the degree of compression depends 154 on the mounting technique used [14]. The initial geometry of the compressed embryo 155 was based on segmented egg shapes obtained from microscopy data [15], while the 156 uncompressed shape was estimated by fitting an ellipsoid with the same volume [41].

157
The spindle offset α for each cell division was statistically sampled from experimental 158 data in compressed conditions from Fickentscher and Weiss [6]. Furthermore, we assume 159 fixed spindle skew angles φ, with values from experiments of Pimpale et al. [39]. In 160 simulations, we assume that cell-cell adhesion is small relative to cortical tension in 161 early C. elegans development [42,43]. A full description of the computational model of 162 C. elegans development is provided in the SI.

163
With these parameters, together with well-characterized cell division timings, we We hypothesized that the effect of confinement alone is sufficient to explain the 175 different volume ratio observed in experiments [15]. To this end, we adopted spindle 176 positions from experiments in confined conditions, but simulated development in an 177 uncompressed egg geometry. However, we observed that simulation results were 178 significantly different from experimental observations, see Fig. 4B. These results suggest 179 that apart from cell shape changes, egg geometry also influences spindle positioning, as 180 the shape changes alone were found to be insufficient to induce the experimentally 181 observed differences in volume. Since spindle positions seem to be the dominant factor 182 in determining daughter volume ratio, our model suggests that any significant change in 183 the volume of daughter cells should be accompanied by a shift in spindle position.  Table 2. This yields estimates of α for the compressed 192 conditions that agree with experimental measurements [6,44], and produces a 193 computational prediction for the adjusted spindle offset in uncompressed conditions. In 194 the latter case, the model predicts a mean absolute shift in α of ±300 nm. The   [15,39]. A linear fitting procedure was used where α values were estimated via interpolation, see S2.
Finally, a sensitivity analysis was performed with respect to the influence of cortex 198 mechanics on cell division asymmetry in the C. elegans embryo, thereby assessing the

208
The effect of relative cortex viscosity κ was found to be negligible, further supporting 209 the notion of increased stability due to geometrical confinement. Similarly, varying the 210 normalized contact friction λ produced no significant effects on cell volume ratio. This 211 latter observation is to be expected, as up to the 6 cell stage embryo, cell motility is low 212 and no drastic cell-cell contact area changes occur.

214
In this work, we presented a computational mechanical model of cytokinesis, which 215 explicitly represents the contractile furrow ring that drives ingression of a cortical shell 216 under mitotic tension. We demonstrated that this model well describes the dynamics of 217 furrow ingression, on the condition that the furrow ring is mechanically dominant 218 compared to the cortex surface tension and effective viscosity. These results agree with 219 the predictions of Turlier et al. who used a visco-active membrane theory of the cortex 220 to provide a general framework for interpreting and characterizing constriction failure 221 and furrow ring closure as a result of underlying furrow ring and cortex mechanics [22]. 222 Furthermore, we extended the cytokinesis model to account for asymmetric cell division 223 by varying the position and angle of the mitotic spindle complex. In case division 224 mechanics are dominant over cortex mechanics, the model reproduces a simple geometric 225 relationship between spindle offset and daughter volume fraction. When confined 226 between parallel plates, daughter cell asymmetry increases when the spindle orientation 227 is parallel to the plate and decreases when it is perpendicular to the plate. For uniaxial 228 compression, these differences in sensitivity can be traced back to a change in the 229 characteristic length in the direction perpendicular to the division plane. For non-trivial 230 compressed states, e.g., cells during C. elegans embryogenesis, similar effects occur 231 where the division direction and characteristic length heavily influence the robustness of 232 the system to any perturbation in spindle positioning. For C. elegans specifically, this 233 partially explains the observed changes in volume ratio for ABa and ABp daughter cells. 234 We demonstrated that that the model is able to recapitulate C. elegans 235 embryogenesis up to the six cell stage, using spindle position and angle as model input. 236 In the case of C. elegans, we conclude that spindle position alone predicts the observed 237 volume fraction after cell division very well. This supports the proposition that during 238 cytokinesis, the mechanical properties of the contractile ring are dominant with respect 239 to cortical properties such as mitotic tension and viscosity. In contrast to our initial mechanical stress, but also provides supplementary stability for asymmetric divisions. 255 We did not investigate the effects of cell-cell adhesion, cell-egg adhesion and 256 embryo-egg volume ratio, as Giammona et al. have already covered these topics [28].

257
Finally, we also did not investigate the important effect of confinement, cell-cell 258 connectivity and signaling on spindle offset α and skew angle φ [10,45]. Incorporating 259 the underlying mechanisms that mechanically regulate these structures will prove 260 crucial for expanding this model towards fully predictive simulations of development.

261
Nonetheless, this work shows the feasibility for the DCM framework to simulate 262 developing systems such as early C. elegans morphogenesis, while allowing for detailed 263 comparison to experimental observations.

265
We approximate the mechanical behavior of the cell cortex as a viscous shell, 266 represented by a triangulated surface mesh where nodal positions act as the relevant 267 degrees of freedom, see Fig. 2A. The cortical shell is parameterized by a thickness t, 268 surface tension γ and viscosity η. Cell volume is conserved using an apparent bulk 269 modulus K . For simulation of multicellular C. elegans development, adhesive cell-cell 270 and cell-egg contact with adhesion energy w (assumed γ) and wet friction constant ξ 271 is implemented [32][33][34]. A full mechanical description of the model can be found in S1. 272 Computational model of the furrow ring 273 Structurally, the furrow ring is composed of fixed length sub-units, in which tension is 274 generated due to depolymerization of sub-units in the presence of end-tracking 275 cross-linkers [17,46]. To model this process, the furrow ring is represented as an elastic 276 ring in the mother cell cortex, consisting of sub-units with a given initial length.

277
Depolymerization is modeled by decreasing the resting length of the elastic sub-units, 278 locally generating tensile stress in the cortex [17]. We do not model furrow-microtubule 279 interactions, but introduce an effective sub-unit depolymerization velocity as with base depolymerization velocity v d = l0 τ d , given initial sub-unit lengths l 0 = ±600 nm and base ring closure time τ d [17,36]. Integrating eq. 3 yields eq. 1, as fit on experimental data by Carvalho et al., when fixing R t = 3.6 µm and R f = 0.8 µm.
Using a linear elastic material model for the furrow ring, the in-plane elastic energy is expressed as ij represent the current and resting distance between nodes 281 ij and where k el is the spring stiffness. The spring force between connected nodes is 282 then given by Consequently, the nodal spring forces are F el i = F el ij and F el j = −F el ij . As nodal 284 positions do not necessarily align with sub-unit endpoints, the furrow ring is 285 approximated by a closed loop consisting of node-node connections closest to the 286 division plane. For a given mesh refinement, this means that the selected node-node 287 edges can be either larger or smaller than l 0 . To address this scaling issue, the effective 288 node-node depolymerization velocity is determined for each unique connection, with scale factor The spring force F el ij is added to the equation of motion presented in S1. Additionally, 291 the local viscosity is increased accordingly, see S1, eq. 5. The local force balance induces 292 March 15, 2022 9/14 cortex deformation, resulting in furrow ingression, see Fig. 1. Cell abscission is triggered 293 when the furrow perimeter reaches a set threshold value, 2πR f . During this phase, 294 furrow ring components are disassembled and the mother cell is split into two new 295 entities, where the daughter cells inherit the mechanical properties of the mother cell.

296
Simulated constriction velocities were estimated based on numerical differentiation of 297 the ring perimeter with a step size of 10 s, while the total cytokinesis duration was 298 estimated by taking the difference between cell abscission and cytokinesis initiation 299 timings.

300
Stabilization of polar instability in dividing cells

301
Asymmetric division of a viscous tensile shell is mechanically inherently unstable due to 302 polar contractility, resulting in the bubble-instability. As we do not account for 303 microtubule-cortex interactions, we use the approach proposed Sedzinski et al., where a 304 linear bulk elastic resistance K div is introduced to account for these interactions, 305 assuring cell shape stabilization and avoiding daughter-half collapse during 306 cytokinesis [31]. This stabilization model assumes cytosolic flow between daughter 307 halves to be tightly regulated so that daughter halves can be treated as individual 308 volumes. Using a proportional pressure controller the reference volume V D of daughter-half D, based on the division plane and mother 310 cell polarization, can be enforced. P div D is integrated over the nodal Voronoi area S i and 311 added to equation of motion, S1, eq. 3. The net reaction force due to the pressure 312 difference between daughter halves, is applied as a reaction force to the furrow ring. This ensures that a proper force 314 balance is maintained at the global (cell) level.