In order to quantify development over the whole tissue from cell to tissue level, we applied our formalism to measure $\mathbf{G}$, $\mathbf{D}$, $\mathbf{R}$, $\mathbf{S}$, $\mathbf{A}$ using 2 h sliding window averages in cell patches of initially 40$\times$40 μm${}^2$, typically encompassing tens of cells that were tracked between 14 and 28 hAPF (Video 4a-e). Dilation and CE rate maps vary smoothly in both space and time, while remaining symmetric relative to the midline axis. This indicates that the averaging length scale (40$\times$40 μm²) and time-scale (2 h) are appropriate to describe the dynamics of the tissue; moreover, the symmetry indicates that each process is regulated. The results were also plotted as 14 h average maps on the last image analyzed to quantify the growth and morphogenesis of each patch between 14 and 28 hAPF (Figure 3f–j and Figure 3—figure supplements 1, 2). We find that tissue dilation mainly occurs through divisions, cell size changes and apoptoses, the dilation rates of the other processes being negligible (Figure 3—figure supplement 1). While the dilation rates are critical to study important processes such as apical constriction or local tissue growth, we focus here on the CE rates to illustrate how the formalism enables to better understand the role of cell processes in morphogenesis, with a focus on oriented cell divisions. The tissue CE rate $\mathbf{G}$ map demonstrates that the tissue undergoes various CE movements, in particular in its posterior medial and lateral domains (Figure 3f, boxes). As expected, the cell division orientation (${\mathbf{D}}_{\mathbf{o}}$) and division CE rate ($\mathbf{D}$) maps (Figure 3e,g) show strong similarities (alignment = 0.94). Importantly, the formalism now enables to compare directly the division CE rate $\mathbf{D}$ to the other CE rates. First, the $\mathbf{G}$ and $\mathbf{D}$ CE rates are roughly aligned in the medial posterior region, while they have clearly different orientations in the lateral domains (Figure 3f,g). Second, both cell rearrangements ($\mathbf{R}$) and cell shape changes ($\mathbf{S}$) have strong CE rates in many domains where cell divisions are also oriented (Figure 3g–i). Third, while cell delaminations are spatially regulated and numerous (~2.5 10³), their CE rates ($\mathbf{A}$) are negligible, and their effect is mostly isotropic (Figure 3j and Figure 3—figure supplements 1e, 2e). Therefore, tissue CE mostly occurs through oriented divisions, rearrangements and cell shape changes. Together this shows how our formalism enables to express quantitatively the relevant information extracted from millions of cell contours over space and time and to separately measure and represent the CE rates of the tissue and of each cell process at the scale of groups of cells in a whole epithelial tissue.

Video 4

Download asset

In order to determine the average and the variability of cell dynamics of a tissue, we developed a general method to register in time and space, and to rescale movies obtained from different animals. This method uses reliable landmarks (for space registration) and prominent, stereotyped rotational movements (for time registration); and it can be applied to mutant conditions and to other tissues (see Materials and methods and below for mutant conditions). Upon space-time registration, each deformation rate measured in each cell patch in 5 hemi-notum movies can be averaged to yield average maps of an archetype hemi-notum representing the dynamics of ~7 10⁶ cell contours, i.e. ~23 10⁶ links, 40 10³ divisions, 5.4 10³ delaminations (Figure 4a,b and Figure 4—figure supplement 1a). For each process and in each region of the archetype hemi-notum, the corresponding biological variability is measured by quantifying its variations between hemi-nota in each region. This was used to determine whether each CE rate was significant in a given region (plotted in color or grey, respectively). Collectively, these approaches yield average maps of the tissue CE rate $\mathbf{G}$ and of the main cell process CE rates ($\mathbf{D}$, $\mathbf{R}$, and $\mathbf{A}$) that can be superimposed for comparison (Figure 4a,b). These maps confirmed that division CE rates are parallel to the tissue deformation rates in some regions but not all (compare regions 1, 2 and 3), suggesting that oriented divisions may not have a single and simple effect in morphogenesis.

Figure 4 with 1 supplement see all

Download asset Open asset

Studying the morphogenesis of epithelial tissues requires the analysis of both the alignment and the magnitude of each cell process CE rate with respect to the local tissue CE rate. This can be achieved by projecting each CE rate ($\mathbf{D}$, $\mathbf{R}$, $\mathbf{S}$ and $\mathbf{A}$) onto the direction of the local tissue CE rate (noted ${\mathbf{u}}_{\mathbf{G}}$) to determine their components along the direction of $\mathbf{G}$, thereby yielding ${\mathbf{D}}_{\mathrm{}//}$, ${\mathbf{R}}_{\mathrm{}//}$, ${\mathbf{S}}_{\mathrm{}//}$ and ${\mathbf{A}}_{\mathrm{}//}$ components respectively (Figure 4c). The projection of $\mathbf{G}$ along its own direction yields the local magnitude of tissue morphogenesis (${\mathbf{G}}_{\mathrm{}//}$). Each of these components along $\mathbf{G}$ can therefore be interpreted as the effective participation of each process in tissue morphogenesis of the wild-type tissue; it should not be confused with a functional role of a cell process that can only be studied using loss and gain of function approaches and modeling. A CE rate aligned with the tissue CE rate has a positive component along tissue morphogenesis, whereas a CE rate displaying a bar orthogonal to the tissue CE rate has a negative component along tissue morphogenesis (Figure 4c). In a tissue where tissue deformation is solely associated with cell divisions, cell rearrangements, cell size and shape changes and apoptoses, $\mathbf{G}$ is equal to the sum of the CE rate of divisions ($\mathbf{D}$), rearrangements ($\mathbf{R}$) and cell shape changes ($\mathbf{S}$), as well as apoptoses ($\mathbf{A}$). The magnitude of local tissue morphogenesis (${\mathbf{G}}_{\mathrm{}//}$) can therefore also be expressed as the sum of the local components of each cell process along $\mathbf{G}$, namely:

Maps of the components associated with each cell process can be then plotted using circles that are hollow or filled for positive or negative components, respectively, and the areas of which scale with the magnitude of each components (Figure 4d–g and Figure 4—figure supplement 1b,c).

We briefly illustrate here how the component measurements can further be analyzed to uncover novel interactions between tissue morphogenesis, cell divisions and other cell processes. The measurements show that the cell division component (Figure 4e) can either be strongly positive in some regions (regions 1 and 2) or strongly negative in others (region 3). In contrast, cell rearrangements and cell shape changes mainly have positive components along tissue morphogenesis (Figure 4f,g). In regions 1 and 2, the positive component of cell divisions corresponds to the one described in the literature, namely that cell divisions and tissue elongation are oriented in the same direction. Conversely, our measurements in region 3 show that cell divisions have a negative component corresponding to divisions taking place mostly orthogonally to the direction of tissue elongation. The measurements of cell rearrangements and cell shape components show that the elongation of the tissue in this region results from the strongly positive components of cell shape changes or of cell rearrangements (Figure 4f,g). Finally in region 4, although most cells divide once, the cell division component is small relative to that of cell shape changes that almost completely accounts for tissue deformation in that region (Figure 4d–g).

The existence of these distinct relationships between process CE rates and total CE rate can be further analyzed by plotting each relative component versus time to characterize its temporal evolution throughout tissue development (Figure 4h–k). In addition to the spatial heterogeneity of morphogenesis, these plots reveal that morphogenesis also strongly varies over time in a given region of the tissue and that it can be further decomposed in different periods (see Figure 4h–k legends for details). More generally, such analyses will provide important insight about the respective roles of each cell process in tissue morphogenesis and their time evolution.

To illustrate the generality of our approach and to determine whether cell divisions can be oriented perpendicularly to the tissue CE rate in another tissue, we performed a similar characterization of the morphogenesis in another epithelium, the pupal wing (Figure 5, Figure 5—figure supplement 1, Video 5a). The important deformations of cell patches over time show that wing morphogenesis is heterogeneous as well, displaying strong tissue deformations in the anterior and posterior domains of the wing, and milder deformations in the medial domain (Figure 5a, Video 5b). To characterize quantitatively these deformations of the tissue and relate them to cellular processes occurring in this wing, we applied our formalism to this pupal wing between 15 and 32 hAPF and determined maps of averaged tissue CE rates $\mathbf{G}$ and cell processes CE rates $\mathbf{D}$, $\mathbf{R}$, $\mathbf{S}$ and $\mathbf{A}$ (Figure 5b–e), as well as their projections onto $\mathbf{G}$ direction, ${\mathbf{G}}_{\mathrm{}//}$, ${\mathbf{D}}_{\mathrm{}//}$, ${\mathbf{R}}_{\mathrm{}//}$, ${\mathbf{S}}_{\mathrm{}//}$ and ${\mathbf{A}}_{\mathbf{}\mathbf{/}\mathbf{/}}$ (Figure 5f–i). Even on a single wing, all CE rate maps display heterogeneous but smooth patterns over time, showing that this scale of space-time averaging (40$\times$40 μm², 2 h) is also suitable for the wing (Video 5c–f). These measurements confirmed the observed heterogeneity in tissue morphogenesis by evidencing major tissue CE rates in anterior and posterior domains of the wing, with a tilt of about ±45° with respect to the horizontal axis, respectively, and smaller CE rates in the medial domain (Figure 5b,f). The CE rates of cell processes $\mathbf{D}$, $\mathbf{R}$, $\mathbf{S}$ display smooth and heterogeneous maps while displaying different patterns in space (Figure 5c–e), while the CE rate of $\mathbf{A}$ is negligible (Figure 5—figure supplement 1b). Like in the dorsal thorax, rearrangements and cell shape changes mainly have positive components along tissue CE rate and they make up most of tissue morphogenesis (Figure 5f,h,i). Like in the dorsal thorax, oriented cell divisions in the wing display more variety in their component along tissue morphogenesis than $\mathbf{R}$ and $\mathbf{S}$: $\mathbf{D}$ has positive components along $\mathbf{G}$ in anterior and posterior wing domains, where tissue morphogenesis is the strongest, while it has negative components in the medial wing domain, where morphogenesis is milder (Figure 5c,g box). Together our analyses of morphogenesis in the notum and the wing illustrate how our formalism enables the quantitative characterization of morphogenesis in various epithelial tissues.

Figure 5 with 1 supplement see all

Download asset Open asset

Video 5

Download asset

ubi-E-Cad:GFP (Oda and Tsukita, 2001) and E-Cad:GFP (Huang et al., 2009) were used for live imaging of apical cell contours in notum and wing pupal tissues. In all experiments the white pupa stage was set to 0 h after pupa formation (hAPF), determined with 1 h precision. For notum tissues, pupae were prepared for live imaging as described in (David et al., 2005). Pupae were imaged for a period of 15–26 h, starting at 11–12 hAPF, with an inverted confocal spinning disk microscope (Nikon) equipped with a CoolSNAP HQ2 camera (Photometrics) and temperature control chamber, using Metamorph 7.5.6.0 (Molecular Devices) with autofocus. Movies were acquired at 25±1°C. We took Z-stacks with 18 to 28 slices (0.5 μm/slice) to make sure we included the adherens junctions. Maximum projections of Z-stacks were obtained using a custom ImageJ routine (publicly available as the ‘Smart Projector' plugin) and have been used for tissue flow analysis. Multiple position movies were stitched using a customized version of the ‘StackInserter' ImageJ plugin. Filming 10 to 12 positions with 0.32 μm spatial resolution (pixel size) every 5 min yielded a tiling of a half dorsal thorax (3 movies), 24 positions yielded a tiling of the whole dorsal thorax (1 large movie, available as Supp. Video in [Bosveld et al., 2012]), resulting in 5 hemi-notum movies.

For trbl^up notum tissues, temporal control of gene function was achieved by using the temperature sensitive Gal4/GAL80^ts (McGuire et al., 2003). Tribbles was overexpressed using an UAS-Trbl transgene (Grosshans and Wieschaus, 2000) during pupal stage. Embryos and larvae were raised at 18°C. Late third instar larvae were switched to 29°C. After 24 to 30 h, pupae were examined. Those which were timed as 11$\pm$1 hAPF were mounted for live imaging at 29°C. Five hemi-notum movies were acquired.

For wing tissue, E-Cad:GFP pupa was prepared for live imaging as described in (Classen et al., 2008). Immersol W 2010 (Zeiss) was used instead of voltalef oil. Pupal wing was imaged for a period of 17 h at 5 min interval, starting at 15 hAPF with an inverted confocal spinning disk microscope (Olympus IX83 combined with Yokogawa CSU-W1) equipped with iXon3 888 EMCCD camera (Andor), an Olympus 60X/NA1.2 SPlanApo water-immersion objective, and temperature control chamber (TOKAI HIT), using IQ 2.9.1 (Andor). Other details of treatment and quantification of wing data will be published elsewhere.

The local tissue flow within notum tissues was first estimated by image cross-correlation velocimetry along sequential images using customized Matlab (Bosveld et al., 2012) routines based on the particle image velocimetry (Rael et al., 2007) toolbox, matpiv (http://folk.uio.no/jks/matpiv). Each image correlation window was a 32$\times$32 pixels square (~10$\times$10 µm², pixel size 0.32 μm), had a 50% overlap with each of its neighbors. The resulting estimate of the velocity field was used to initiate the tracking procedure, and also for time registration of different movies.

Images were denoised with the Safir software (Boulanger et al., 2010). Cell contours were determined and individual cells were identified using a standard watershed algorithm. Errors were corrected through several iterations between manual and automatic rounds. Automatic rounds consisted of a custom automatic software (Matlab) which tracked segmented cells through all images and detected each event where two cells fuse (Bosveld et al., 2012). The tracking relied on the comparison between cells in one image and the following, based on overlap of cells as well as centroid positions. Divisions were detected and we checked the sizes and elongation of daughter cells. Delaminations were characterized by cell size decrease before cell disappearance, and distinguished from fusions. A fusion between times $t$ and $t+\delta t$ is an artifact which almost always reveals an error of segmentation which is either a false positive cell junction at time $t$, or a cell junction missing at time $t+\delta t$. Both by manual tests on small subsamples, and by automatic tracking of false cell appearances or disappearances, we estimated the final relative error rate to be below ${10}^{-4}$, which was sufficient for the analyses presented here. The whole notum movie contains ~8.8 10³ cells at the beginning. After one to three divisions per cell, several delaminations, and a flow of cells out of the field of view, the final image contains ~18.4 10³ cells. Altogether, on the five hemi-notum movies, ~7.7 10⁶ cell contours were segmented and tracked. For trbl^up tissues, there were ~3.7 10⁶ cell contours. Cells were larger and easier to segment, resulting in errors even smaller than in wild-type tissues (data not shown). Cells in the anteriormost part get out of focus due to tissue flow and elongation, so that we cannot track them throughout the whole movie; similarly, on lateral sides, some cells become visible during the movie (purple cells in Video 3). We observed a total of ~40 10³ divisions and ~5 10³ delaminations in the wild-type tissues, ~1.7 10³ divisions and ~7 10² delaminations in trbl^up tissues. This yielded maps of cell apical surface area and shape, cell centroid displacements (Videos 1 and 3), cell lineages, and neighbour changes.

For wing tissue, ~2 10⁶ cell contours were segmented and tracked. The segmentation was less manually corrected than in the notum. Errors in segmentation or tracking resulted in cell fusions, cell integrations and fluxes; their measured values were small enough that the measurement for other cell processes were not significantly affected. Interestingly, this shows the robustness of our formalism with respect to errors in input data. Altogether, ~13 10⁶ cell contours (~40 10⁶ cell-cell junctions) were analyzed.

In epithelial tissues where cells are confluent, assuming mechanical equilibrium, cell shapes are then determined by cell junction tensions $\gamma$ and cell pressures $p$, and information on force balance can be inferred from image observation (Brodland et al., 2010; 2014; Ishihara and Sugimura, 2012; Chiou et al., 2012; Ishihara et al., 2013). For instance, if three cell junctions which have the same tension end at a common meeting point, their respective angles should be equal by symmetry, and thus be 120° each. Reciprocally, any observed deviation from 120° yields a determination of their ratios. The connectivity of the junction network adds redundancy to the system of equations, since the same cell junction tension plays a role at both ends of the junction.

Mathematically speaking, there is a set of linear equations, which involve all cell junction tensions, and cell pressure differences across junctions. We simultaneously estimate tensions and pressures by using Bayesian statistics formulated by maximizing marginal likelihood, or equivalently, by minimizing the Akaike Bayesian information criterion (ABIC) (Akaike, 1980; Kaipio and Somersalo, 2006). Our expectation that junction tensions are distributed around a positive value is incorporated as a prior (Ishihara and Sugimura, 2012; Ishihara et al., 2013). This inverse problem requires to invert large matrices, whose typical size is a few ${10}^4\times {10}^4$. Our custom program takes advantage that these matrices are sparse (http://faculty.cse.tamu.edu/davis/suitesparse.html), which not only increases the speed of resolution, but also minimizes the ABIC more robustly. It infers all junction tensions up to only one unknown constant which is the tension scale factor, and an additional unknown constant which is the average cell pressure. By integrating tensions and pressures, their contributions to tissue stress can be calculated in any given cell patch $\mathcal{\wp }$, again up to these two constants, with the Batchelor formula (Batchelor, 1970; Ishihara and Sugimura, 2012)

where $\mathbb{𝟙}$ is the two-dimensional identity matrix, $\mathcal{𝒜}$ is the cell patch area, ${\mathcal{𝒜}}_c$ the area of cell $c$, the sum is over each cell $c$ in the patch and its neighbours $k$, ${\overrightarrow{ℒ}}_{ck}$ is the vector representing the chord of the junction between cells $c$ and $k$ (it links two vertices, its orientation being unimportant). Only the junction tensions contribute to the deviator of $\mathbf{\sigma }$ (hence called ‘junctional stress' in the main text), and thus determine the main direction and the difference of two eigenvalues of $\mathbf{\sigma }$.

Regarding this junctional stress, we recall three points:

• Junctions contribute to transmit stress between cells even if the biological origin of this stress lies in non-junctional cytoskeletal elements and other structures in the body of the cell. Forces created in the cell body or at the cell boundary are transmitted to its neighbours by contact at the cell-cell interface. These forces transmitted by contact can be decomposed in components parallel and perpendicular to the junction, namely tension and pressure, which we both estimate. Changes of cell volume lead to measurable changes in pressure and thus to detectable changes in our stress estimation.

• We focus on junctional stress because of the growing consensus that it is a important determinant of tissue morphogenesis, as measured by laser ablations of (and outside of) adherens junctions (for review see [Lecuit et al., 2011]), or directly probed by experimental manipulation of adherens junctions and compared with simulations (Bambardekar et al., 2015). Note that we cannot measure other contributions to the total stress, for instance through cell membrane parts far from adherens junctions.

• The force inference method used to estimate the junctional stress is independent of other contributions to stress or of external force contributions; it thus remains valid even in the presence of cell-substrate interactions that can cause an additional contribution to the total mechanical interactions.

Each of the five wild-type hemi-notum movies was oriented with the midline along the top horizontal side. The whole notum movie was cut into two hemi-notum movies; the left one was flipped to be oriented like the others movies. As spatial landmarks, we use the macrochaetae, which we identify manually at ∼16h30 APF $\pm$ 10 min, when the sensory organ precursor (SOP) cells divide (Gho et al., 1999), and the following results were independent of the choice of this reference time. For each movie, each macrochaeta is labelled $p$, with $p=1$ for the closest to the midline, and the largest $p=7$ or 8 according to the movie. For each of the five movies, labelled $a=1$ to 5, we call $x$ axis the midline oriented from posterior to anterior, $y$ axis the perpendicular axis oriented from medial to lateral, and as origin $O$ of $(x,y)$ axes the barycenter of the five macrochaetae closest to the midline, $p=1$ to 5. A box is defined on the first image as the set of cells which barycenter is within a 128$\times$128 pixels square (~40$\times$40 μm²). A box contains several tens of cells, and has a 50% overlap with each of its neighbors. Each box is then labelled by two integer numbers $(m,n)$ which define our two-dimensional space coordinates on a fixed square lattice (Eulerian representation). Averaging over space is implemented by averaging over $(m,n)$. When the whole notum movie is analyzed alone as a whole, the midline is chosen as $x$ axis.

Data near the sample boundary are less reliable because of poorer statistics. In order to improve the quality of all quantitative calculations and visual representations, a cell patch at location $(m,n)$ at time $t$ near the boundary of the sample was assigned a weight $0\le w^a(m,n,t)\le 1$, as follows.

We defined the ‘bulk’ of the tissue in Eulerian description as boxes containing only core cells, and in Lagrangian description as patches placed at least three patches away from the boundary.

In Eulerian description, we defined the ‘relative area’ $a_r(m,n,t)$ of a box as the sum of the surface area actually occupied by recognizable cells in this box, divided by the surface area of the box. It was close to 1 in all boxes of the bulk, possibly slightly larger than one due to large cells with their barycenter in the box but some of their area outside of the box, and possibly lower than one due to junction pixels (junctions are one pixel wide, they belong to the box but are not assigned to a cell). In Lagrangian description, $a_r(m,n,t)$ was the number of cells in the box, divided by its maximum value in the tissue: in the first images of the movie it was slightly less than 1 over the whole bulk, while it decreased by a factor at least 2 towards the end of the movie.

Its minimum value $a_{\mathrm{bulk}}\left(t\right)$ in the bulk at time $t$ was used to normalized all values as

which can occur only for boxes out of the bulk, and

which occur for all boxes in the bulk. This way, $a_{\mathrm{norm}} (m,n,t)$ is equal to 1 all over the bulk, and gradually vanishes when approaching the tissue boundary as cell patches contain fewer cells.

Since the results presented below relied on the inversion of the texture matrix (Section C.1.1), we measured the condition number of the texture matrix in this box. The condition number of a function with respect to an argument measures how much the output value of the function can change for a small change in the input argument. For a symmetric 2$\times$2 matrix such as $\mathbf{M}$, the condition number is $\parallel \mathbf{M}\parallel \parallel {\mathbf{M}}^{-1}\parallel$, always larger than 1. We used its inverse, Matlab’s ‘reciprocal condition number’ $Rc$, which is 1 when $\mathbf{M}$ is isotropic for instance and vanishes when $\mathbf{M}$ is not invertible (i.e. for a box which contains only one half-link, or half-links which are all parallel to each others). Its minimum value $R{c}_{\mathrm{bulk}}\left(t\right)$ in the bulk at time $t$, typically $\sim$ 0.5, was used to normalized all values as

This way, $R{c}_{\mathrm{norm}}(m,n,t)$ is equal to 1 all over the bulk, and vanishes when approaching the tissue boundary. The weight of the box was the square of the normalized relative area times normalized reciprocal condition number:

Note that since the stress estimation does not use the inversion of the texture matrix, the weights used for stress measurements are simply $\left[a_{\mathrm{norm}}\right(m,n,t){]}^2$.

In this work, we mostly consider quantities that have been averaged over some time period ($\Delta t=$2h or 14h) rather than their instantaneous values that can be noisy. For any quantity at $(m,n,t)$, calling $q(m,n,t)$ its instantaneous value, its weighted average over time is $Q(m,n,t)$ and is defined as follows:

the sum being made between $\tau =t-\Delta t/2$ and $\tau =t+\Delta t/2$. The corresponding mean weights averaged over the same period are:

with $N_{\Delta t}=\Delta t/\delta t$ is the number of frames during $\Delta t$, the time between two frames being $\delta t$. Weights $W^a(m,n,t)$ decrease according to the number of data in the box, from 1 for a box fully in the bulk of the tissue between $t\pm \Delta t/2$, to 0 for a box outside of the tissue boundary between $t\pm \Delta t/2$. These weights were used in all calculations, averages and graphs. It was even used in maps as an opacity, so that data become more transparent near the boundary where they are less reliable (Figure 3, Figure 3—figure supplements 1, 2, Figure 5, Figure 7 and Videos 4, 5).

We define the scalar product of two second-order tensors $\mathbf{Q},{\mathbf{Q}}^{\text{'}}$ in dimension $d$ as follows:

The scalar product of $\mathbf{Q}$ with itself is the square of its norm

Note that with this definition, the identity $1$ is unitary:

In dimension 2, for two deviatoric tensors represented as bars, $\mathbf{Q}.\mathbf{Q}\text{'}=\parallel \mathbf{Q}\parallel \parallel \mathbf{Q}\parallel \text{'}\cos \left(2\mathrm{\Delta }\theta \right)$ where $\mathrm{\Delta }\theta$ is the difference in their bar directions. $\mathbf{Q}.\mathbf{Q}\text{'}$ is positive if the bar eigendirections are aligned (close to 0°), negative if these directions are perpendicular (close to 90°), and vanishes in between (close to 45°). The alignement coefficient between two tensors is:

It is close to 1 if both tensors are aligned everywhere, close to $-1$ if if both tensors are perpendicular everywhere, and vanishes if the tensor directions are independent. Equation 8 is an average of $\cos \left(2\mathrm{\Delta }\theta \right)$ with weights combining $w^a(m,n)$ and the bar lengths.

We translated the five movies in order to match their origins $O$. The position of macrochaetae is reproducible from one animal to the other: their standard deviation is 5.5 μm along $x$ and 6.3 μm along $y$. We define as an archetype of species $s$ (e.g. the wild-type) a virtual animal which macrochaeta $p$ would be the barycenter of the actual macrochaetae $p$ in the five movies:

We then rescale each actual movie $a$, separately along $x$ and $y$ axes. This is necessary at least for movies taken at different resolutions, or for mutants. The procedure is as follows. Along $x$, the multiplicative factor which minimizes the dispersion of macrochaetae of $a$ from the archetype is:

In the wild-type tissues, we find ${\alpha }^a$ ranging from 0.96 to 1.07, with a standard deviation of 0.03. Along $y$ axis, we perform the same analysis, and also include the position of the midline as an independent information; we find rescaling factors of 1 $\pm$ 0.04. We thus observe that for these five hemi-notum movies the rescaling is not crucial. The variability in tissue scale for these wild-type tissue movies is lower than the variability in macrochaetae positions (not shown). After rescaling, the residual dispersion is slightly lower than the initial one: the standard deviation of macrochaetae position becomes 4.5 μm along $x$ and 5.5 μm along $y$. The referential $(O,x,y)$ in the archetype defines the grid, whereby (0,0) is centered around one box. In the trbl^up tissues, standard deviations of macrochaetae positions were 6.2 μm along $x$ and 6.3 μm along $y$. With respect to the archetype, the standard deviations were 10.2 and 7.9 μm, respectively. With the same procedure as the wild-type tissue, they were rescaled; after rescaling, the standard deviations were 6.2 and 5.7 μm, respectively. Note that the change in temperature from 18°C to 29°C has apparently no effect on the tissue shape, according to tests performed on the posterior part of the wild-type tissue (Bosveld et al., 2012).

While the hAPF was determined with 1 h absolute precision, the tissue rotation rate analysis provided a better relative precision that we used to synchronize the different movies, as described in detail in (Bosveld et al., 2012). In the region of the tissue located near the origin of axes, we observed that the rotation rate systematically passed through a maximum during the development: this rotation peak could be used as a biological reference time. For that purpose, the rotation rate was measured and spatially integrated over a rectangular reference window. Plotting this average versus frame number yielded a bell-shaped curve (shown in Figure S4 of [Bosveld et al., 2012]). We applied a time translation to each movie so that these curves overlapped. We matched the portion of the curve which had the steepest slope: we thus used the time corresponding to 3/4 of the peak value, in the ascent (rather than the maximum itself, which by definition had a vanishing slope). Its average value was 18:40 hAPF. Hereafter, ‘hAPF' indicates the time after this temporal translation has been applied. For instance, after synchronization, the maximum of the contraction-elongation and rotation rates were consistently found at 19:20 and 19:40 hAPF, respectively. This determination reached a $\pm$ 1 interframe (i.e. $\pm$5 min) relative precision. After this synchronization, we determined that the period of time included in all movies was 13:55-27:55 hAPF (hereafter and in the main text, rounded to ‘14-28 hAPF’), corresponding to 169 frames, 68 interframes analyzed in what follows. Global time averages were performed over all these frames. Sliding window averages were performed over 2 h every 5 min. The macrochaetae were again manually determined at the time 16h30 APF determined with this synchronization. To improve the precision, the time and space registration was iterated: it changed only slightly, evidencing the robustness of this double registration. The trbl^up tissues were synchronized with the same procedure as the wild-type tissue, up to a $\sim$10 min precision. For experiments performed at 29°C, the development is accelerated. This was taken into account by dividing time intervals by 0.9 (as determined with $\pm 0.05$ precision by the widening of the wild-type tissue rotation peak [Bosveld et al., 2012]).

These spatial and temporal adjustments allowed us to assign a system of space-time coordinates $(m,n,t)$ common to the 5 hemi-notum movies, from pupae with a given species $s$ (wild-type or mutant). We again checked, now with spatial details, the reproducibility from one pupa to another. To estimate the averages and standard deviation among animals at each space-time point $(m,n,t)$, for each movie labelled $a=1$ to 5, we used the time-average weights $W^a(m,n,t)$ defined in Equation 4. For any given time-average quantity $Q^a(m,n,t)$, we defined its ensemble average ${\overline{Q}}^s(m,n,t)$ over the $N_s$ movies of the species (here, the wild-type genotype, $N_s=5$), local in space and time, as

The corresponding mean weights averaged over the same set of animals are:

These weights were used in maps as an opacity so that measurements become more transparent near the boundary where data come from fewer animals and are less reliable (see Figure 4, Figure 6, Figure 6—figure supplement 1, Figure 7).

The biological variability of the quantity $Q(m,n,t)$ was determined by measuring the weighted standard deviation $\delta Q(m,n,t)$ for the species

The fraction in Equation 12 generalizes the usual term for unbiased standard deviation, ${\left(N_s-1\right)}^{-1}$, which is recovered when all weights $W^a$ are equal to 1. This unbiased weighted standard deviation can be calculated for each of the different independent components of a symmetric tensor $\mathbf{Q}$. In general, these components are $\delta (Q_{xx}^s+Q_{yy}^s)$, $\delta (Q_{xx}^s-Q_{yy}^s)$ and $\delta Q_{xy}^s$. When we consider only the deviator of $\mathbf{Q}$, the components are $\delta (Q_{xx}^s-Q_{yy}^s)$ and $\delta Q_{xy}^s$. We compare them to the biological variability $\delta Q$, and define that a measurement at a given position and time $(m,n,t)$ is significant if it satisfies:

namely if:

Using the significance criterion of Equation 14, measurements larger than the biological variability are plotted in color while smaller ones are shown in grey. Note that if the direction of the anisotropic part of the tensor varies a lot between animals, this variation decreases the amplitude of the anisotropic part of the mean tensor $\overline{Q}$, or equivalently the amplitude of $\overline{Q}$, and therefore the measurement is considered as non significant based on Equation 14.

We numerically simulated the different processes using the cellular Potts model (Graner and Glazier, 1992; Glazier and Graner, 1993), an algorithm relevant in biology to describe variable cell shape, size, packing and irregular fluctuating interfaces of cells in 2 or 3 dimensions (Marée et al., 2007; Bardet et al., 2013). Each cell is defined as a set of pixels, here on a 2D square lattice; their number defines cell area. The pixelisation of the calculation lattice can be chosen to match the resolution of experimental images. A cell shape changes when one of its pixels is attributed to another cell instead. We used periodic boundary conditions, with external medium surrounding cells (a state without adhesion or area and perimeter constraints). We use here a simple version where cells minimise their surface energy. Energy minimisation uses Monte Carlo sampling and the Metropolis algorithm, as follows. We randomly draw (without replacement) a lattice pixel and one of its eight neighboring pixels. If both pixels belong to different cells, we try to copy the state of the neighboring pixel to the first one. If the copying decreases the total energy, we accept it, and if it increases the total energy, we accept it with a probability <1 known as ‘fluctuation allowance’.

Divisions were implemented as follows. Cells were growing with an initial asynchrony in their cycle ranging from 0 to 40% (to avoid divisions all occurring at the same time). As they grew, the group of cells underwent a force gradient along the $x$ axis which tended to stretch each individual cell shape as well as the cell patch shape along the $x$ axis. Cells divided only once when their target area had doubled, along their long axis; they did not regrow after division, thereby recovering their initial size. Cells divided along their long axis, which resulted in divisions oriented along $x$. Cells tended to relax their elongated shapes before or after divisions, resulting in some rearrangements. Each image was 800$\times$800 pixels² and was cropped to yield the final movie. There were 50 Monte-Carlo steps between two successive images.

Other processes were implemented in a similar way (except when stated otherwise), as follows. Stretching the pattern and allowing for cell shape relaxation led to oriented rearrangements. Affinely stretching the pattern by direct image manipulation using an image treatment software instead of Potts simulations (hence not followed by any cell shape relaxation) led to strong cell shape changes. Delaminations were obtained by gradually decreasing the cell target area. Cell integration movie was produced by reversing the order of images of delamination simulation movie. Fusions were forced on cells by random removal of cell-cell junctions (the same image was gradually Gaussian-blurred then thresholded), and cells were then let to relax their shapes. Boundary flux was implemented by gradual removal of successive external cell layers. Rotation of the cell patch was achieved by direct image manipulation. A simulation movie is made with 241 frames, 240 interframes. It would be analogous to a experimental movie of 240 $\times$ 5 min, or 20 h. Hence an experimental scale of 10^-2 h^-1 would be equivalent to 8.3 10^-4 interframe^-1 in simulations.

The approach which leads to the formalism and to the decomposition into separate cellular processes is described in details in the Appendix. We describe here its practical implementation.

In practice, when studying a monolayer of cells, we are mainly interested in morphogenetic movements within the plane of the monolayer: we focus on in-plane components and actually implement the formalism in two dimensions, as follows. We use two successive images of the segmented and tracked movie. In both images, for each cell $c$ the list of neighbour cells $k$ is identified. A half-link $ck$ is the link between a cell $c$ center and its neighbour $k$ center, and is listed independently from the half-link $kc$ oriented from $k$ to $c$. In both images, we list all half-links. The rare cases where four cells meet (along a vertex, in 2D, or along a line, in 3D) are listed separately: in what follows, it is possible to treat them separately if needed. This is what we do, and at the end we lump their contribution with that of other cells.

The tracking enables to classify each half-link into one and only one category (which ensures the completeness of the formalism): appeared or disappeared through one of the processes, or conserved between both images. For instance, when a cell divides, all its half-link disappearances and the appearances of half-links for its daughters are included in the division process (conversely, the ‘division orientation' includes only the link between the newly created daughter cells) (Figure 1a,b). Similarly, when cell $c$ undergoes a delamination, some of its former neighbors enter in contact: their new half-links are also counted in the delamination contribution (Figure 1a,b). The half-link $kc$ can be in a different category from half-link $ck$, for instance if in this interframe cell $c$ undergoes a division while cell $k$ undergoes a delamination, rearranges with other neighbours, or exits at the sample boundary. The result is a classification of all half-links in each image.

All individual movies were analyzed using Lagrangian measurements, thanks to the tracking of cells and patches illustrated in Figure 3d, Figure 5a, Videos 1 and 3. To coarse-grain the measurements, the whole image is subdivided into cell patches $\mathcal{\wp }$. In the Lagrangian description chosen here, each patch is tracked from one frame to the next. After a division, daughter cells are assigned to the patch of their mother. All cells and links belonging to each patch are used to calculate tissue deformation and the deformations associated with each cell processes occurring between two images, and we then turn them into rates and average them over time (see Appendix, sections B and C). Those deformation rates are related by the balance equation (see Equation 77):

We measure separately each term of both sides of Equation 15, by assigning each link geometric change to $\mathbf{G}$, and each link topological change to a specific cell process $\mathbf{P}$, taking advantage of their completeness. Finally, we systematically check that Equation 15 is satisfied.

We consider any of the cell process (e.g. divisions, rearrangements, cell shape and size changes, delaminations, … ), which we note $P$, measured by the tensor $\mathbf{P}$ (e.g. $\mathbf{D}$, $\mathbf{R}$, $\mathbf{S}$, $\mathbf{A}$, … ). While the isotropic part of $\mathbf{P}$ is a scalar, which directly quantifies the growth rate, the CE part of $\mathbf{P}$ is a tensor characterized by an amplitude and a direction that differ from the amplitude and direction of tissue CE rate $\mathbf{G}$. In order to determine the effective contribution of process $P$ to local tissue deformation, we project it onto $\mathbf{G}$. To do so, in each patch from the tissue, we first define the unitary tensor ${\mathbf{u}}_{\mathbf{G}}$ that is aligned with $\mathbf{G}$ (Figure 4c):

which trivially ensures $\parallel {\mathbf{u}}_{\mathbf{G}}\parallel =1$. Then we project process $P$ contribution $\mathbf{P}$ along ${\mathbf{u}}_{\mathbf{G}}$ using the scalar product defined in Equation 5 (Figure 4c):

We call this projection ${\mathbf{P}}_{//}$ the component of $\mathbf{P}$ parallel to tissue morphogenesis (or ‘along $\mathbf{G}$', for short). It is expressed as a rate of change per unit time, i.e. hour^-1. It represents the effective contribution of $\mathbf{P}$ to tissue morphogenesis. The additivity (Equation 15) also applies separately to these components:

In 2D, ${\mathbf{P}}_{//}=\parallel \mathbf{P}\parallel \cos \left(2\Delta \theta \right)$, where $\Delta \theta$ is the difference in directions of $\mathbf{G}$ and $\mathbf{P}$ bars. Thus, if a process CE rate $\mathbf{P}$ is exactly parallel to the tissue CE rate $\mathbf{G}$, it has a positive component on $\mathbf{G}$, which is exactly its own amplitude: ${\mathbf{P}}_{//}=\parallel \mathbf{P}\parallel$. If the $\mathbf{P}$ bar is rather parallel to $\mathbf{G}$, it has a positive component ${\mathbf{P}}_{//}$ on $\mathbf{G}$. If the $\mathbf{P}$ bar is at 45° angle to $\mathbf{G}$, its component ${\mathbf{P}}_{//}$ on $\mathbf{G}$ is zero. If the $\mathbf{P}$ bar is rather perpendicular to $\mathbf{G}$, its CE rate has a negative component ${\mathbf{P}}_{//}$ on $\mathbf{G}$. If the $\mathbf{P}$ bar is exactly perpendicular to $\mathbf{G}$, its CE rate has a negative component on $\mathbf{G}$, which is exactly minus its own amplitude: ${\mathbf{P}}_{//}=-\parallel \mathbf{P}\parallel$. As a consequence, the component of tissue CE rate along itself, ${\mathbf{G}}_{//}$, is the amplitude of $\mathbf{G}$, and is positive by construction. Figure 4c shows the sign of ${\mathbf{P}}_{//}$.

In this section, we list the various sources of uncertainties that we encounter, present our methods to decrease them as much as possible, and discuss why their impact on our analyses remains limited. They fall into two main categories: those related to the acquisition of data which will be used as input for the formalism; and those related to the formalism itself.

• Errors due to image analysis affect the input of our formalism. Segmentation errors can result in false identification of cells, while tracking errors can result in false identification of cells and their lineages. It is thus necessary to choose an image acquisition time sufficient to ensure an image contrast and quality which enable a segmentation with a low error rate. This sets a minimal value to the acceptable time difference $\delta t$ between two successive images. In the notum and wing tissues we observe here, fusion and integration events are used as markers of segmentation and/or tracking errors. In the notum, we use them as feedback to correct the segmentation until the contributions of fusions and integrations are negligible with respect to the contributions we want to measure.

• The tissue is a three-dimensional (‘cuboïdal’) monolayer and we do not assume it to be 2D. However, we perform its study in 2D. More precisely, we image it using the E-Cadherin:GFP marker, which labels the apical adherens junctions. Moreover, several quantites we measure are actually 2D, such as the velocity field, or its gradient which represents the morphogenesis. Other quantities like cell division orientation or stress field are 3D but can be studied in 2D independently from what occurs in the third dimension. Finally, some 3D processes like delamination in which apical area shrinks have an effect on tissue morphogenesis which is similar to that of an apoptosis where a cell entirely shrinks, and is thus not necessary to distinguish within the scope of the present study.

• The 3D structure of the notum also plays a role because the tissue is curved. It is imaged in 3D using a confocal spinning disc microscope; the height of its surface is recorded and we can entirely reconstruct its curvature. We can then obtain the local angle $\theta$ of the tissue surface with the horizontal plane: in the notum, this is mainly along the direction $y$ (medio-lateral), since along the axis $x$ the tissue curvature is much smaller. We obtain the projection factor cos$\theta$ which should in principle be applied as a correction in the corresponding direction either to the raw image before any treatment, or a posteriori to results of the formalism. We measure $\theta$ is at most 0.01 radian over the tissue and reaches at most 0.02 radian in the most lateral part of the whole thorax image. The absolute correction factor is the same for all tensors, and it is negligible with respect to 1 in the results presented here (1-cos$\theta$ ~10^-4), so that it does not affect significantly the results presented here. For simplicity, we choose not to perform any correction to our results, knowing that we can perform the correction a posteriori if it appears necessary in the future. All these effects are completely negligible in the wing, which is much flatter than the notum.

• For a given set of data used as input to our formalism, the formalism itself does have uncertainties which affect its output. For instance, the time difference $\delta t$ between two successive images should be small enough so that the fraction of conserved links remains larger or comparable to the fraction of non-conserved links, and that our measurement of tissue growth and morphogenesis (which is based on conserved links) accurately quantifies the actual morphogenesis. Since $\delta t$ should be sufficiently large to keep a good image quality, the question is how to choose an optimal $\delta t$. We find that in our case, $\delta t=5$ min is optimal, and for these value both constraints are satisfied: the image is good enough to be automatically segmented for a large part, and the morphogenesis is well quantified by our measurement. Finally, in order to validate our workflow, we have checked that with the data of the preceding paper (Bosveld et al., 2012) our new formalism recovers consistent results. In conclusion, our pipeline which links the image acquisition to the data representation (Figure 1e) is valid and its accumulated uncertainties are lower than the biological variability within an animal or between different animals.

In mutant tissues, the space and time registration are performed as in the wild-type tissues. Mutant tissues can exhibit a different variability than wild-type tissues; in practice, we observe a larger variability in trbl^up mutant tissues than in wild-type tissues. Moreover, wild-type and mutant tissues can be registered together. The total morphogenesis as well as the measurement of each cell-level process is computed similarly in a mutant tissue. We compare them term by term with the corresponding values in wild-type tissues. We assay at which time and position the mutation has a significant effect with an inter-genotype variability larger than the intra-genotype one. It is also possible to subtract term by term each measurements performed in the wild-type tissues ${\mathbf{P}}_{\text{wt}}$ and the measurements performed in the mutant tissues ${\mathbf{P}}_{\text{mutant}}$. We can plot the measurement difference $\mathrm{\Delta }\mathbf{P}$ defined as $\mathrm{\Delta }\mathbf{P}={\mathbf{P}}_{\text{mutant}}-{\mathbf{P}}_{\text{wt}}$ this difference represents the part of ${\mathbf{P}}_{\text{mutant}}$ that has been added by the mutant condition to ${\mathbf{P}}_{\text{wt}}$. It represents the effect of the overexpressed gene in the process $P$ of wt morphogenesis, and its projection $\mathrm{\Delta }{\mathbf{P}}_{//}$ along the wild-type tissue CE rate represents its effective contribution to wild-type morphogenesis.

Once the cells have been segmented and tracked throughout the movie (see ‘Materials and methods’ for details), our goal is to extract useful quantitative information on divisions, rearrangements, cell size and shape changes and apoptoses, and quantify the relative importance of these cell-scale events in complex morphogenetic processes. Several descriptive terms, such as cell or tissue polarity, oriented cell divisions, convergence-extension, or intercalation, have an underlying common point: their quantitative description refers to an amplitude, an anisotropy and a direction. The mathematically robust objects that encompass simultaneously these informations are the tensors, which also offer the advantage to be valid in spaces of two, three or any dimensions. Using tensors is particularly important in a heterogeneous tissue like the notum, where morphogenetic flows of various amplitude, direction and anisotropy occur at different locations.

For instance, counting the divisions results in integer numbers, and statistics on cell division directions results in angle distributions. While these classical descriptions of the cell division rate and orientation are essential, we need to quantitatively determine to which extent they contribute to tissue morphogenesis. Characterizing divisions using tensors can result in a better description of oriented cell divisions and a better understanding of their effects. Similarly, if a cell rearrangement is followed immediately by its inverse (which often occurs either in reality, or as an artifact when segmenting a short junction), this leads back to the initial state, and their opposite effects cancel out; while if two rearrangements occur in the same direction, their effects cumulate. Counting cell rearrangements as integer numbers leads to $1+1=2$ rearrangements in both examples, while a tensor analysis better captures in each case its actual impact to morphogenesis.

There is a need for a formalism with the following requirements : it should include a self-consistent, rigorous definition of the tensors used for measurements; enable a decomposition of tissue deformation into the sum of contributions coming from elementary cell processes; the tensors should be significantly measurable with a good signal-to-noise ratio, and visually representable; the formalism should be purely descriptive, in the sense that it should quantify observed morphogenetic changes independently of the biological or physical origin of the forces that drive these changes, whether they are internal forces such as active cytoskeletal forces, viscous dissipation or cytoplasm pressure, or external forces such as interaction with another cell layer or solid substrate.

We aim here to generalize and take further the approaches used in four articles which have begun to apply this type of methods to developmental biology, enabling quantitative comparison between experiments and simulations:

• Rauzi et al. (2008) used a matrix called ‘texture’ as defined in Equation 6 of (Graner et al., 2008) to quantify cell shape changes (see their Figures 5a and S4). It is based on the links connecting each cell center with the centers of its neighbors: it expresses, in μm², the variance of link length in each direction. Separately, rearrangements were counted and quantified as integer numbers.

• Aigouy et al. (2010) used a former version of the texture as defined in the Equation 1 of (Aubouy et al., 2003), which expressed the variance of the cell junctions in each direction. More precisely, they called ‘cell elongation’ (defined in their Eq. S25) the deviatoric part of the texture. Rearrangements were quantified separately, and as integer numbers (see their Figure 4).

• Butler et al. (2009) used their own method (introduced in [Blanchard et al., 2009]). They computed the contribution of cell shape changes to flow. They then subtracted the cell shape change rate from the tissue shape change rate. They thus indirectly inferred the contribution of all types of topological changes to morphogenesis as tensors, i.e. including their amplitude and their direction.

• In a preceding paper (Bosveld et al., 2012), we adapted the formalism of Graner et al., 2008, already validated on non-biological cellular networks (Cheddadi et al., 2011), to quantify tensorially the contributions of cell shape changes and rearrangements to tissue shape changes in the Drosophila scutellum, a subpart of the dorsal thorax. We took advantage of averages over several animals to improve the signal-to-noise ratio, and to subtract a mutant from a wild-type tissue, evidencing the effect of the mutation and determining at which place the difference was significant. In (Bardet et al., 2013) we used a similar approach in an intervein of the Drosophila wing.

Two additional formalisms have also been proposed recently. Economou et al. (2013) use scalars (rather than tensors) to analyze the one-dimensional elongation of a developing mouse tissue. Etournay et al. (2015) is based on the estimation of the local tissue deformation at the sub-cellular scale of triangles formed by three contiguous links surrounding each tricellular junction (three-fold vertex) and analyses the pupal wing development in details, using a tensorial approach.

The formalism presented here is tensorial, to take into account the variable orientations of cell processes. It accounts for processes like cell divisions, apoptoses or delaminations which can drastically modify cell number in a proliferative tissue. We made it complete by including the remaining processes that can occur in other tissues, or because of segmentation mistakes, or simply because of the limited field of acquisition of the microscope: coalescences (fusions) and new cell integrations, as well as fluxes of cells through the boundary of the tissue of interest or of the microscope field of view.

In a cellular material such as a living tissue, the material states can be characterized by the links joining the centers of mass of neighboring cells, thereby tiling the tissue with an irregular triangular networks of links (Figure 1, Figure 1—figure supplement 1, Video 1). Then, by keeping track of those links over time, the geometrical and topological changes in the tissue can also be characterized, and related to the elementary cell processes (cell kinematics) (Video 1) (Graner et al., 2008; Etournay et al., 2015; Tlili et al., 2015).

Therefore, as a starting point, we use a statistical definition similar to (Graner et al., 2008) to characterize the state and the deformation of the material network of links. Each change regarding a link (appearance, disappearance, change of orientation or length) is assigned to one and only one process, guaranteeing that the morphogenetic events are decomposed into individual processes and that these processes add up to tissue growth and morphogenesis. One of our main goals is to obtain a balance equation relating the local tissue deformation to the elementary cell processes making up this deformation, thereby connecting the cell and tissue scales.

To quantify tissue deformation and all cell processes, we aim at defining dimensionless tensors in order unify their characterizations, but also to relate their statistical descriptions to quantities used in continuum mechanics. Indeed, the continuum mechanics description facilitates comparison between different experiments, or between experiments, simulations and theory. It requires to treat the tissue as continuous, that is, describe cell patches without any explicit signature of each individual cell size. Cell processes are described statistically rather than with details. For each cell process, it is possible to construct a biologically relevant continuous counterpart of the statistical quantity, namely a quantity which is dimensionless (or a rate), with no size.

The procedure presented here improves and generalizes the ones presented in (Graner et al., 2008; Bosveld et al., 2012; Bardet et al., 2013). One of the main difficulties to overcome is that biological processes such as divisions or delaminations can change drastically the number of cells and links on the image in proliferative tissues. If we had directly adapted the formalism presented in the above papers (Section A.2), the tensor describing cell shape and size changes would have been strongly affected by the changes in cell and link numbers, and would thus have lost a clear interpretation. In addition, divisions and delaminations would have been mixed with rearrangements.

The formalism presented here provides a unified quantitative characterization of tissue deformation and of all cell processes making up this deformation. All these quantities are characterized with dimensionless symmetric tensors that quantify the tissue strain and the elementary strains associated with each cell process, independantly of rigid body movements. It provides a simple balance equation where the tissue strain is equal to the sum of the cell process strains. Each tensor remains biologically interpretable as a separate contribution to the tissue morphogenesis. Importantly, the formalism is valid at arbitrarily large and heterogeneous deformations, whether the tissue is highly proliferative or not, and it is defined in a space of arbitrary dimension, and is therefore directly implementable in two or three dimensions.

After having briefly recalled some basic definitions of deformation in continuum mechanics and of the statistical description of a cellular material that constitute our starting points, we define a coarse-grained deformation gradient tensor $\mathbf{F}$ for a general deformable cellular material (section B). Then we define a coarse-grained strain tensor to quantify the total strain in the cellular material ${\mathbf{G}}^{*}$, and we relate it to the strains associated with elementary cell processes (${\mathbf{S}}^{*}, {\mathbf{P}}^{*}$) with a simple and general balance equation (section C). Incidentally, this enables us to define the material geometrical (${\mathbf{S}}^{*}$) and topological (${\mathbf{T}}^{*}$) strains. Finally, we compare our approach with the ones previously published (section D).

In a cellular material, the material state can be characterized by the discrete links joining the centers of mass of neighbor cells. Together, these links constitute a network of irregular triangles (Video 1c). By following those links over time, the material deformation can be characterized as well, and related with elementary cell processes (Video 1c–h) (Graner et al., 2008; Etournay et al., 2015; Tlili et al., 2015).

During growth and morphogenesis, the tissue patches change in shape and size, and their deformation can be characterized in situ by the changes in length and direction of the links that are conserved between successive images.

We now show how we can relate continuum mechanics quantities and the statistical description of a deformable cellular material by building a deformation gradient tensor $\mathbf{F}$, not at the cell scale, but rather directly coarse-grained at the cell patch scale (based on the changes of the material links between two states) to quantify the deformation of the material, and by relating it to $\mathbf{G}$.

In this section, we relate the newly defined tissue strain tensor ${\mathbf{G}}^{*}$ (Equation 63) to all elementary cell processes, namely cell size and shape changes, divisions, rearrangements, apoptoses/delaminations, integrations, fusions, and inward/outward fluxes. This will enable the definition of meaningful tensors to characterize quantitatively the strains associated with each cell process.

First, now that we have shown that one can build a proper strain tensor to characterize tissue deformation based on $\mathbf{G}$, we rewrite Equation 34 using Equation 32, and we use Equation 53 to make all symmetric tensors dimensionless:

Using Equation 63 to make ${\mathbf{G}}^{*}$ appear, we get:

Now we have a first balance equation linking the tissue strain ${\mathbf{G}}^{*}$ to tensors quantifying cell processes. Indeed, $\mathrm{\Delta }\tilde{\mathbf{M}}$ is related to the cell size and shape changes in the patch between the initial and final states (Graner et al., 2008), while the tensors $\tilde{\mathbf{P}}$ quantify the changes in the link network topology due to each cell process $P$. However, to have tensors that meaningfully and unambiguously characterize the contribution of each elementary cell process to tissue strain requires some additional steps detailed below.

When the number of half-links in the patch varies from $N$ to $n$ between the initial and final states, the variation in texture has two entangled contributions, one due to the cell size and shape changes, and one due to the changes in number of half-links. They can be disentangled by rewriting $\mathrm{\Delta }\tilde{\mathbf{M}}$ as follows:

The first term in Equation 67, by renormalizing the final texture $\tilde{\mathbf{m}}$ by the ratio $\frac{N}{n}$, quantifies the cell size and shape changes independently of the changes in half-link number in the patch. The second term in Equation 67 involves the total difference in half-link numbers in the cell patch between the initial and final states, $\mathrm{\Delta }N=n-N$. Importantly, this difference can be decomposed into the sum of the changes in half-link number due to each elementary cell process $P$ that impacts the network topology $(\mathrm{\Delta }{N}_{\mathrm{P}}=n_{\mathrm{P}}-N_{\mathrm{P}})$:

$\mathrm{\Delta }{N}_{\mathrm{P}}$ therefore quantifies the contribution of each topological process to the total change in half-link number within the cell patch. $\mathrm{\Delta }{N}_{\mathrm{P}}>0$ for divisions, integrations or inward fluxes since those cell processes correspond to a net creation of half-links; $\mathrm{\Delta }{N}_{\mathrm{P}}<0$ for apoptoses/delaminations, fusions or outward fluxes since those cell processes correspond to a net loss of half-links; and $\mathrm{\Delta }{N}_{\mathrm{P}} \approx 0$ for rearrangements since they mostly conserve the number of links, the only exception involving cells at the image boundaries.

Combining Equations 67,68 one can therefore rewrite Equation 66 as:

where each $\frac{\mathrm{\Delta }{N}_{\mathrm{P}}}{n}\tilde{\mathbf{m}}$ term related to the change in half-link number due to each topological process $P$ has been gathered with its corresponding tensor $\tilde{\mathbf{P}}$ in the sum, thereby gathering all terms impacted by topological changes.

We can now define our final expressions for the tensors representing the strains of elementary cell processes contributing to tissue strain:

• the first term of Equation 69 defines the tensor quantifying the contribution of changes of cell size and shape to the tissue strain ${\mathbf{G}}^{*}$. It represents the geometrical strain in the tissue and is called ${\mathbf{S}}^{*}$ (see section C.3.1).

• the second term of Equation 69 represents the total topological strain in the tissue and is called ${\mathbf{T}}^{*}$. Each term of the sum over topological processes defines the tensor ${\mathbf{P}}^{*}$ quantifying the contribution of the cell process $P$ to total tissue topological strain ${\mathbf{T}}^{*}$ (see section C.3.2).

With those new notations Equations 32,69 therefore read:

The total tissue strain tensor ${\mathbf{G}}^{*}$ can therefore be expressed as the sum of the geometrical strain tensor ${\mathbf{S}}^{*}$ and of the total topological strain tensor ${\mathbf{T}}^{*}$. The latter can be further broken down into elementary topological strains tensors ${\mathbf{P}}^{*}$ associated with each cell process changing the link network topology. We detail and comment these new definitions in the next sections.