Active cell divisions generate fourfold orientationally ordered phase in living tissue

Morphogenesis, the process through which genes generate form, establishes tissue-scale order as a template for constructing the complex shapes of the body plan. The extensive growth required to build these ordered substrates is fuelled by cell proliferation, which, naively, should destroy order. Understanding how active morphogenetic mechanisms couple cellular and mechanical processes to generate order–rather than annihilate it–remains an outstanding question in animal development. We show that cell divisions are the primary drivers of tissue flow, leading to a fourfold orientationally ordered phase. Waves of anisotropic cell proliferation propagate across the embryo with precise patterning. Defects introduced into the nascent lattice by cell divisions are moved out of the tissue bulk towards the boundary by subsequent divisions. Specific cell proliferation rates and orientations enable cell divisions to organize rather than fluidize the tissue. We observe this using live imaging and tissue cartography to analyse the dynamics of fourfold tissue ordering in the trunk segmental ectoderm of the crustacean Parhyale hawaiensis beginning 72 h after egg lay. The result is a robust, active mechanism for generating global orientational order in a non-equilibrium system that sets the stage for the subsequent development of shape and form.


Initial lattice construction
For the division wave/random division simulations, the initial lattices were generated by first constructing a perfect hexagonal lattice with a specified length and height.Row identities were extracted from this initial stage and a subset of the rows in the middle of the lattice were selected to be a part of the division order generated using the row identities.For the random division simulations, this initial division order was then randomly permuted.All simulations were run with an initial set of five dividing rows.Gaussian white noise with a signal-to-noise ratio of 18 was then added to cell positions to produce a nonuniform hexatically ordered initial state (i.e.quasi-long range sixfold orientational order and no translational order).The decision to use a hexatically ordered initial state was due to the necessity of having an initial configuration with well defined rows and with no fourfold orientational order.Properties of the initial lattices are displayed in Fig. S19.After the initial lattice was constructed, the parameter A 0 was set to be the mean area of all the cells in the lattice and P 0 = 4A 0 .Lattices were constructed so that all cells had an approximate initial length scale = 1.

Growth, division, relaxation, and confinement
After the initial lattice was constructed, simulations begin with the first division wave.The cell that is next to divide in the division order undergoes a series of growth steps.In each growth step, the parameters P 0,c and A 0,c were increased, fixing P 0,c = 4A 0,c .Lattice connectivity was then updated by implementing T1 transitions for all edges below a user Algorithm 1 Vertex Model Simulation Algorithm Require: Initial cell lattice, initial division order for first wave Relax initial cell lattice for Each division wave do Update division order for current wave for Each cell in the division order do for Number of growth steps do Update A 0 , P 0 for the next cell in the division order Update mesh connectivity due to T1 transitions Relax cell lattice Enforce confinement end for Re-set A 0 , P 0 for the cell about to divide Update mesh connectivity due to T1 transitions Randomly choose division axis orientation Divide cell along division axis Relax cell lattice Enforce confinement end for end for specified length threshold.Boundary edges whose shared external angle dropped below a user specified threshold were then fused together to prevent self-interpenetration of the tissue and other physically erroneous configurations.For all of our simulations, we set the T1 edge length threshold min = 0.05, the edge merger threshold θ min = 20 • , and the number of growth steps N growth = 50.The geometric parameters P 0,c and A 0,c were increased so that A 0,c = 1.5 A 0 after N growth growth steps.Additionally, for all of our numerical experiments T1 transitions were prohibited in the bulk.
After all parameter updates and topological corrections, the lattice was then relaxed to an equilibrium configuration.The energy was minimized with the GNU Scientific Library Multidimensional Minimizer, using the Polak-Ribiere conjugate gradient algorithm.The wrapper to this minimization procedure was extended from the code used in [1].Following relaxation, confinement was enforced by clipping the x-coordinates of lattice vertices to remain within a user specified domain, equal to the width of the original lattice domain plus one cell length.This choice was consistent with the observed behavior of the germband during the first stage of growth (Fig. 4C).Using a large number of growth steps ensured that no more complicated constraint enforcement procedure was necessary.
After all of the growth steps were concluded, the cell was then forced to divide.Division orientations were drawn from a circular von Mises distribution given by where µ is the center of the distribution, k is its concentration, and I 0 (k) is a modified Bessel function of the first kind.We set µ = π/2 for all of our simulations, i.e. a vertically oriented division along the A-P axis.See Fig. 5F for a visual depiction of the distributions for all values of k used in the simulations.Following division, the tissue was relaxed and confined again and the whole procedure continued for the next cell in the division order.Once all of the cells in the division order had divided, the division order was updated for the subsequent wave so that the new order would also proceed as a division wave.We only simulated a single division wave in our numerical experiments.Visualization and order parameter calculations were performed on the Voronoi tessellation of the centroids of each cell polygon in order to be consistent with the methodology used to analyze the data.

Fluid mechanical model
We investigated the role of cell divisions in producing tissue scale flow fields through the use of an active hydrodynamic model.Continuum mechanics is an effective method for capturing how force balance is associated with local deformation and flow [2].In the simplest models, the forcing that generates flow can be neatly decomposed into a set of separate passive and active contributions.The tendency of cells to resist deformation and the viscous coupling between cells translate local forces into long range flows.Interactions among these various contributions, both locally and over large scales, can be then be built up with increasing degrees of complexity to create more accurate descriptions of biological realities.This framework has proved effective in numerous biological applications elucidating how active forces generate non-trivial flow fields [3][4][5].
The tissue scale mechanics of epithelia is best understood as emerging from the collective mechanics of the tissue's constituent cells [6].The cell scale mechanics are dominated by the influence of intracellular cytoskeletal cortices [7,8].These structures are capable of supporting stresses throughout the cell's bulk interior and can deform the cell either actively or in response to external stresses.The individual cortices of neighboring cells are coupled together by cadherin mediated adherens junctions [9] into a global trans-cellular mechanical network.Crucially, these networks are both 'active' and 'adaptive'.Active contributions, such as cell divisions [10] or the contraction of trans-cellular actomyosin networks [4], act both locally and over large scales to deform the tissue.These active deformations generate passive stress distally from the activity via the coupling between cells.Stresses created in this fashion can then be relaxed either locally, via internal rearrangement of the cytoskeletal cortex [11,12], or at the tissue scale, via cell rearrangements [3].The result is an exotic type of viscoelasticity, wherein the tissue responds elastically to active stresses over short timescales and then reverts to a fluidlike flow at longer timescales as the tissue adaptively relaxes stress.When the timescale of growth (i.e.rate of cell division, etc.) is long compared to the timescale of mechanical relaxation, the tissue essentially behaves like a slowly creeping fluid in quasistatic mechanical equilibrium with both the active and external forces.For our purposes, a detailed description of the complex cellular processes mediating the active stresses and adaptive relaxation is unimportant.Instead, we adopt an approximate continuum scale description of the tissue that captures the relevant behaviors of the system.
At short time scales, the tissue behaves elastically [13,14].Incremental increases in strain generate corresponding increases in stress.These stresses are then subsequently relaxed over a time scale τ R .The total stress σ = σ e + σ a can be decomposed into a sum of the elastic stresses σ e and the active stresses σ a .For simplicity, we assume that the tissue is a flat two-dimensional surface.Supposing the tissue behaves like a Maxwell viscoelastic fluid, the time-evolution of the elastic stress is given by where σe = ∂ t σ e + v • ∇σ e denotes a convective time derivative, v is the local velocity or rate of displacement, and µ and λ are the Lamé parameters characterizing a linear isotropic stress-strain relationship.The Lamé parameters are assumed to be spatially homogeneous, but may depend on time.A more complete description would include the corotational time derivatives of σ e on the left hand side, but we omit them for simplicity.The first two terms on the right hand side describe the generation of stress in proportion to the rate of strain.The final term parameterizes the relaxation of the elastic stress.We have made the assumption that σe only depends on σ a implicitly through the final term τ −1 R σ e (see the force balance condition below).This is not meant to imply that adaptive cell behaviors do not actually depend explicitly on the active stresses.It is merely a mathematical simplification.In reality, there may exist a complex interplay between mechanics, gene expression/cell fate, and adaptive cellular behaviors, but this lies beyond the scope of this work.
In general, the tissue flow velocity can be described by the Cauchy momentum equation [2] ρ where ρ denotes the material density and F EXT denotes any external body forces acting on the system.If the rate of mechanical relaxation is sufficiently fast to reach quasiequilibrium, σe vanishes, leaving us with Substituting this into Eq.( 5), we find that In other words, the transient elasticity of the system gives rise to fluidlike behavior in the quasi-stationary regime parameterized by a set of effective viscosities.The effective shear viscosity ν 1 ≡ τ R µ tends to resist shearing motion, while the effective bulk viscosity ν 2 ≡ τ R (µ+λ) resists the isotropic compressible part of the flow.If the effective viscosities are sufficiently large and the corresponding motion is sufficiently slow we can neglect the inertial term relative to the viscous terms (ρ v << ν∇ 2 v).If, for simplicity, we furthermore assume that there are no external forces acting on the system, we are left with the following force balance condition defining the tissue flow velocities in terms of the active stresses which can be solved given suitable boundary conditions.

Cell divisions and active stresses
In general, there may be numerous different types of active stresses acting within the tissue during growth [4,10,15].These various sources of active stress will likely not contribute with equal importance to the shaping of the nascent tissue.Our observations of the cell division choreography indicate that cell proliferation is the primary factor mediating tissue velocities.In order to simplify our model, we assume that the only relevant active stresses are due to these cell divisions.Our first consideration in modeling active stresses due to cell divisions is to understand how division events might be incorporated into a spatiotempoally coarse grained scheme.With respect to timing, daughter cell separation in real cell division events occurs over a short, but finite time [16].We ignore the complexity of the rapid subcellular properties mediating mitosis and model cell divisions as instantaneous events (for our purposes any process that occurs faster than 5 min, the time resolution of our microscope data, is considered instantaneous).Real cell divisions push on the surrounding tissue, which responds elastically over the short time scale of the actual division, before remodeling to relax the resulting stress.In our model, both the plastic strain due to a division and the subsequent relaxation are assumed to occur instantaneously so that the fluid is always in quasistatic mechanical equilibrium.With respect to spatial coarse graining, we note that the length scale of a single cell division is small compared to the scale of the entire germband.Numerical experimentation showed that approximating division events as point force dipoles produced inaccurate flow fields.We therefore model divisions as small, but finite size events to regularize this inconsistency.We also make the approximation that the inclusion lives in an infinite, otherwise quiescent medium in order to make analytic progress.
The mathematical structure of the fluid mechanical equations of motion in Eq. ( 8) is identical to the structure of the Navier equations of linear elasticity.Exploiting this similarity, we can model the instantaneous displacement (read velocity) due to a cell division as that of an Eshelby inclusion.In the context of classical elasticity, an Eshelby inclusion is a finite subregion of an elastic body which undergoes a permanent plastic deformation [17].Eshelby inclusions, and the associated theories of elastic multipoles [18], have proved useful in numerous applications including understanding shear localization in amporphous solids [19] and the properties of mechanical metamaterials [20].In the current context of viscoelastic tissue growth, this displacement corresponds to both the plastic deformation and the subsequent viscous relaxation, which are assumed to both occur instantaneously on time scales relevant to growth.Relying on the well-trodden mathematical history of Eshelby inclusions, we may, in the following derivation, use terminology similar to Eshelby's original exposition for elastic materials (e.g.eigenstrain).For clarity, we emphasize again that we are modeling the tissue as a viscoelastic fluid, not an elastic solid, and that the parameters of our model are effective viscosities, not elastic moduli which would vanish in an orientationally ordered phase.
In particular, we choose to model the division as a circular Eshelby inclusion of radius a and eigenstrain * αβ = (M/2) δ αβ + q (2n α nβ − δ αβ ), where δ αβ is the Kronecker delta and n is a unit vector along the division axis.Note that we have adopted an index notation to better account for the high-rank tensorial nature of the following calculations.Greek indices vary in the set {x, y}.We also adopt the Einstein summation notation so that all repeated indices are summed.It was shown by Eshelby that the constraining medium generates a constant strain within the inclusion c αβ = S αβγδ * γδ , where S αβγδ is a constant tensor for any elliptical inclusion.For a circular inclusion which yields the following constrained strain within the inclusion Thanks to the linearity of the system, the corresponding velocity, i.e. instantaneous displacement, within the inclusion can now be computed as The velocity outside of the inclusion will satisfy the biharmonic equation ∇ 2 ∇ 2 v out α (x) = 0, subject to continuity at the surface of the inclusion and must also tend to zero as r → ∞ where r ≡ ||x||.Recalling the radial solutions of the biharmonic equation in 2D (i.e., 1, ln r, r 2 , and r 2 ln r), we construct the general solution by considering all combinations of derivatives of the radial solutions that are linear in the eigenstrain, tend to zero at infinity, and transform like a vector field: where A, B, C, and D are constants to be determined.Recall that the velocity outside of the inclusion is a solution to the equation The fact that ∇ 2 ∇ 2 v out α (x) = 0 is a necessary, but insufficient condition on any solution of Eq. (13).With this fact in mind, we begin to calculate the constants by first computing In the first equality, we used the identity to cancel the terms proportional to A, B, and C. In the second equality, we made use of the fact that ∂ 2 (r 2 ln r) ∂x β ∂x β = 4 (ln r + 1).
Next we compute where we have once again made use of both Eq. ( 15) and Eq. ( 16).Substituting these results into Eq.( 13) yields which, when simplified, reveals that The velocity outside the inclusion can therefore be written as The following identities are now required to make further progress: Substituting these identities directly into Eq.( 20) and simplifying allows us to re-write the velocity in the following form The remaining constants can be computed by enforcing continuity of the velocity field at the boundary of the inclusion where r = a.First, we know that the third term that is cubic in x must vanish when r = a since the velocity within the inclusion is linear in x.
The dependence of v out (x) on C can now be removed yielding (26) Referencing Eq. ( 11), we continue to enforce continuity at the boundary of the inclusion and match the coefficients of the velocities inside and outside the inclusion term by term.The term proportional to * x β tells us and the term proportional to * ββ x α tells us Finally we arrive at the following form for the velocity outside of the inclusion In summary, the instantaneous velocity induced by a cell division is has the following vectorial form or, in terms of the effective viscosities This equation can be non-dimensionalized by introducing the dimensionless parameter in terms of which the division velocity is We note that the parameter ν is not the same as the Poisson ratio introduced in Eq. 11 of the Methods  In general, small differences in orientation will result in significantly lower mean order parameters compared to the magnitude of the single cell order parameters.This effect is especially noticeable in the high order case.).(H) The magnitude of the average fourfold orientational order parameter for each parasegment shifted in time so that T = 0 corresponds to the first cell division of mitotic wave 2. (I) The magnitude of the average sixfold orientational order parameter for each parasegment shifted in time so that T = 0 corresponds to the first cell division of mitotic wave 2.     The orientation of the average fourfold and sixfold order parameters relative to the A-P axis.(E) The Frobenius norm of the average full cumulative strain tensor, the anistropic part of the the average strain tensor, and the isotropic part of the average strain tensor.(F) The Frobenius norm of the full average strain-rate tensor, the anisotropic part of the average strain-rate tensor, and the isotropic part of the average strain-rate tensor.

Fig. S1 .
Fig. S1.Surface extraction via level sets.An illustration of the morphological snakes method used to extract the surface of interest in the growing embryo.(A) The spherical initial condition (zero completed iterations).(B) The segmented surface.

Fig. S4 .
Fig. S4.Effects of averaging the discrete fourfold orientational order parameter.(A and B) Average fourfold order parameter over the third order natural neighborhood of the single cells shown in (Fig. 2, B and C).The dotted orange boundary denotes the seed cell.The solid orange boundary denotes the averaging region.Arrows indicate the orientation of the single cell order parameter for each cell in the averaging region.(A) Average order parameter in a disordered region.(B) Average order parameter in a highly fourfold ordered region.In general, small differences in orientation will result in significantly lower mean order parameters compared to the magnitude of the single cell order parameters.This effect is especially noticeable in the high order case.

3 D
Fig. S5.Trunk ectodermal germband exhibits neither sixfold nor translational/smectic order.Measured values extracted from a single embryo.(A) The two-point sixfold orientational correlation functions at three representative time points.All time points exhibit exponential decay.Vertical dotted line indicates largest lateral length scale in the system.(B) The isotropic pair correlation function.Vertical dotted lines indicate length scale of the system along the A-P and D-V axes.All time points exhibit exponential decay, implying no isotropic translational order.(C -D) Anisotropic variations of the pair correlation function.All time points for both variations decay exponentially.The lack of anisotropic translational order along any preferred direction preclude the possibility of smectic order in the system.(C) The pair correlation function measured along the D-V axis.(D) The pair correlation function measured along the A-P axis.

3 CFig
Fig. S6.Trunk ectodermal germband exhibits neither sixfold nor translational/smectic order (cont.).Orientational order correlation functions and isotropic pair correlation function for a separate embryo at the same representative times.(A) The two-point fourfold orientational correlation functions at two representative time points.Tissue exhibits quasilong range order at the at T = 85.8 h.(B) The two-point sixfold orientational correlation functions at three representative time points.All time points exhibit exponential decay.(C) The isotropic pair correlation function.All time points exhibit exponential decay, implying no isotropic translational order.
Fig. S9.Finite size effects on orientational and translational order.Finite size effects were investigated using synthetic data sets describing both fourfold ordered (A-D) and totally disordered systems (E-H).Systems were generated with square aspect ratios and system sizes L = 10, 20, 40, 80 cell lengths.N = 25 synthetic data sets for each system size.(A) Distributions of the mean global fourfold order parameter | ψ 4 | for an orientationally ordered system.Distributions are much more closely packed around the mean for larger systems, but the mean of the distributions are constant for all system sizes.Circle markers denote distribution medians and bold horizontal lines denote distribution means.(B) Fourfold orientational correlation function for an ordered system.Vertical dotted lines indicate system size.Correlations approach and maintain a constant value over the entirety of the system.(C) Radially symmetric pair correlation function g(r) − 1. Pair correlation function decays algebraically over most of the system size.(D) Horizontal cut of the pair correlation function g(r, x) − 1. Correlations also decay algebraically, albeit extremely slowly.(E) Distributions of the mean global fourfold order parameter | ψ 4 | for a disordered system.The mean value for each distribution is more variable across system size than in the ordered case, but all values are still consistent with a disordered system.Mean values appear to monotonically approach zero as system size in increased.Circle markers denote distribution medians and bold horizontal lines denote distribution means.(F) Fourfold orientational correlation functions decay exponentially for all system sizes.Some erroneous spikes are observed at distances larger than the system size.(G) Radially symmetric pair correlation function g(r) − 1 decays exponentially for all system sizes.(H) Horizontal cut of the pair correlation function g(r, x) − 1 decays exponentially for all system sizes.

Fig. S10 .
Fig.S10.Temporal profile of orientational order is similar across different parasegments.All parasegments belong to the same single embryo.(A to G) Per-parasegment profile of of the magnitude of the average fourfold and sixfold orientational order parameters over time.Ordering of panels (A to G) mirrors anterior to posterior ordering of the physical segments (i.e. the parasegment shown in (A) is anterior to the parasegment shown in (B), etc.).(H) The magnitude of the average fourfold orientational order parameter for each parasegment shifted in time so that T = 0 corresponds to the first cell division of mitotic wave 2. (I) The magnitude of the average sixfold orientational order parameter for each parasegment shifted in time so that T = 0 corresponds to the first cell division of mitotic wave 2.

Fig. S11 .
Fig. S11.Division wave pattern is preserved across embryos.(A to C) The location of mitotic wave division events over time in a different embryo than the one shown in (Fig 3, D to F). Shapes indicate the parasegment within which a division occurs.Indicated lines are linear fits to all division events associated to a particular mitotic wave.(A) shows the location of each division along the A-P axis.The speed of mitotic wave 1 is 8.7 ± 0.8 µm/hr.The speed of mitotic wave 2 (AB) is 7.7 ± 1.6 µm/hr.The speed of mitotic wave 2 (CD) is 8.0 ± 1.2 µm/hr.(B) and (C) show the location of each division along the D-V axis.Division times in (C) have been normalized to the occurrence of the first division event associated with a particular wave in a specific parasegment.The speed of mitotic wave 1 is 12.6 ± 1.1 µm/hr.The speed of mitotic wave 2 (AB) is 10.7 ± 0.5 µm/hr.The speed of mitotic wave 2 (CD) is 10.3 ± 0.6 µm/hr.

Fig. S12 .
Fig. S12.Distribution of cell cycle durations throughout growth.Solid lines are Weibull distributions fit to the cell cycle duration histograms.Average duration and standard deviation for mitotic wave 1 is 10.7 ± 2.2 hours from N = 117 cell cycles.Average duration and standard deviation for mitotic wave 2 is 8.7 ± 1.5 hours from N = 27 cell cycles.Average duration and standard deviation for differential cleavage is 4.6 ± 2.0 hours from N = 7 cell cycles.All cell cycles drawn from the same single embryo.

Fig. S14 .
Fig. S14.The velocity of the mean cell division event.Analysis performed on 190 division events drawn from a single embryo.(A) The orientation of the velocity of the mean cell division event.(B) The x-component of the velocity of the mean cell division event.(C) The y-component of the velocity of the mean cell division event.

Fig. S15 .
Fig. S15.Gradients of the velocity of the mean cell division event.Analysis performed on 190 division events drawn from a single embryo.(A) The divergence of the mean cell division velocity.(B) The 'curl' of the mean cell division velocity.(C) The x-component of the Laplacian of the mean cell division velocity.(D) The y-component of the Laplacian of the mean cell division velocity.

Fig. S20 .
Fig. S20.Orientations of cell divisions are not predicted by local mechanical or geometric signals.Analysis performed on 222 mitotic wave divisions and 67 differential cleavage divisions drawn from a single embryo.(A) Histogram of orientations of cell divisions relative to the principal axis of the strain-rate of the corresponding parent cell integrated over the hour preceding cell division.(B) Histogram of the orientations of cell divisions relative to the axis of elongation of the corresponding parent cell averaged over the hour preceding cell division.(C) Histogram of the orientations of cell divisions relative of the orientation of the fourfold order parameter of the corresponding parent cell averaged over the hour preceding cell division.(D) Histogram of the orientations of cell divisions relative to the orientation of the sixfold order parameter of the corresponding parent cell averaged over the hour preceding cell division.

Fig. S22 .
Fig. S22.Timing of cell divisions are not predicted by local geometric signals (cont.).Analysis performed on 184 cell cycles drawn from a single embryo.(A to F) Various geometric fields averaged over a normalized cell cycle.Separate averages are included for each mitotic wave and differential cleavage.(A) The magnitude of the mean fourfold orientational order parmeter.(B) The orientation of the mean fourfold orientational order parmeter.(C) The magnitude of the mean sixfold orientational order parameter.(D) The orietnation of the mean sixfold orientational order parameter.(E) The mean eccentricity of an ellipse fit to each cell.(F) The mean cell density.

Fig. S23 .
Fig. S23.Timing of cell divisions are not predicted by local mechanical signals (cont.).Analysis performed on 184 cell cycles drawn from a single embryo.(A to H) Various mechanical fields averaged over a normalized cell cycle.Separate averages are included for each mitotic wave and differential cleavage.(A) The Frobenius norm of the mean strain-rate tensor.(B) The Frobenius norm of the isotropic part of the mean strain-rate tensor.(C) The Frobenius of the anisotropic part of the mean strain-rate tensor.(D) The trace of the mean strain-rate tensor.(E) The Frobenius norm of the mean cumulative strain tensor.(F) The Frobenius norm of the isotropic part of the mean cumulative strain tensor.(G) The Frobenius norm of the anisotropic part of the mean cumulative strain tensor.(H) The trace of the mean cumulative strain tensor.