## Abstract

Computer-assisted tracking of the shapes of many cells over long periods of development has driven the exploration of novel ways to quantify the contributions of different cell behaviours to morphogenesis. A handful of similar methods have now been published that are used to calculate tissue deformations (strain rates) in epithelia. These methods are further used to quantify strain rates attributable to each of the cell behaviours in the tissue, such as cell shape change, cell rearrangement and cell division, that together sum to the tissue strain rates. In this review, aimed at developmental biologists, I will introduce the general approach, characterize differences in current approaches and highlight extensions of these methods that remain to be fully explored. The methods will make a major contribution to the emerging field of tissue mechanics. Precisely quantified strain rates are an essential first step towards exploring constitutive equations relating stress to strain via tissue mechanical properties.

This article is part of the themed issue ‘Systems morphodynamics: understanding the development of tissue hardware’.

## 1. Introduction

The drive for new ways to quantify changes in the shapes of tissues in embryo and organ development has been pushed by recent technological advances and pulled by a renewed interest in tissue mechanics. The push results from the avalanche of four-dimensional imaging data and the extraction, using computer-assisted cell-tracking methods, of the trajectories, shapes and connectivity of many cells in two or three dimensions over time [1–6]. These datasets need to be interpreted to extract biologically relevant measures. The first and easiest metric to extract from these datasets is cell velocity, but velocities are generally only of interest in their relative sense, compared to the velocities of other cells or relative to a substrate or landmark. For example, the migration of mesodermal cells of the extending chick axis was initially thought to be directional, leading to a search for guidance cues. Subtracting the velocity of the underlying ECM, the substrate relative to which these cells were migrating, revealed that they were performing undirected random walks, essentially diffusing [7]. For epithelia, raw velocities are probably even less useful, since epithelial cells are constrained in their movement to the layer to which they belong. A quantification of the relative movements of cells should, therefore, be an essential first analysis for epithelial data, and methods to achieve this will be the focus of this review.

The pull, driving the direction in which new methods are focused, is the re-emergence of interest in tissue mechanics [8–11]. Gene expression affects tissue morphogenesis through active force generation, mediated by tissue material properties and by stress imposed by neighbouring tissues and embryo architecture (figure 1*a*). We know about gene expression in often remarkable detail [15,16], while quantitative tissue mechanics had been in relative abeyance since D'Arcy Thompson [17], who connected diverse body shapes through simple geometric transformations, implying differential growth rates. We have begun to explore how genetics and mechanics interact, in order to understand normal and dysregulated development [18–20]. It is not implausible to extrapolate from recent work that certainly some, maybe much, of what we do not understand about developmental dysregulation and cancer could be due to errors in mechanical control. A salutary example of what we have been missing is the ‘infection’ of tumorigenicity in somatic cells by nearby tumours through sustained mechanical pressure, leading to dysregulation of β-catenin signalling [21]. Current efforts to measure epithelial mechanics can be classified into four types: (i) measured by physical perturbation [14,22]; (ii) inferred qualitatively from strain rate patterns [23]; (iii) inferred quantitatively from cell geometries [24–28]; (iv) inferred quantitatively from the dynamics of cell behaviour [13,29]. Though often not explicitly, all approaches assume an underlying tissue rheology that is often expressed as a constitutive relationship between stress and strain, through tissue mechanical properties (e.g. figure 1*b,c*; and see Glossary of terms in box 1). The simplest parameter to measure in any constitutive relationship is the *strain* rate, and this is the most obvious first step in any investigation into tissue mechanics. Unhelpfully, the word ‘strain’ in normal speech describes force. In physics, strain is a deformation resulting from a force. In biology, strain maps are the empirical description of how wild-type and mutant phenotypes arise, through cells and patches of tissue moving during morphogenesis. The question this review addresses is how can we quantitatively resolve complex spatio-temporal strain maps into biologically meaningful cell behaviours?

### Glossary.

*Kinematics* describes the deformations of objects, as distinct from *dynamics*, which analyses the forces that cause them. Both are branches of classical *mechanics*.

A *strain* measures the deformation of an object relative to a reference configuration. A *strain rate* is a change in strain over time.

A *constitutive equation* formalizes the relationship between stress and strain in a material. *Rheology* concerns the relationship between stress and strain in fluid-like matter.

A *tensor* describes a linear geometric transformation, independent of any particular coordinate system. The *trace* of a matrix is the sum of the diagonal elements. For a matrix tensor the trace defines the rate of area change.

A *T1 process* is a plastic change of *topology*, or connectivity, where two cells break planar contact, replaced by a new contact between a different pair. A *T2 process* occurs when a cell leaves or joins a planar array of cells.

*Pure shear* describes an equally balanced convergence–extension strain motif. *Simple shear* describes a gradient of strain perpendicular to the direction of movement, and is equivalent to a pure shear with rotation.

For the most part, this review will concentrate on tissue *kinematics*, the geometric transformations that quantify morphogenetic change without reference to force. I will focus on discrete and continuous methods to measure strain rates in two dimensions (2D), particularly relevant to epithelial morphogenesis and will conclude by briefly revisiting tissue mechanics in the light of good kinematic data.

## 2. Kinematic strain rates

Tissue and component cell behavioural strain rates have been quantified in various mono-layered epithelia of the developing *Drosophila*; the germ-band [1,23,30] and amnioserosa [1,13,31–33], the wing disc [29], the salivary gland placode [34] and the pupal notum [35]. They have also been characterized in the more complicated pseudo-stratified epithelia of the zebrafish trunk ectoderm [1] and early forebrain [36], the developing tooth of the mouse [37] and during the formation of the primitive streak in the gastrulating mouse embryo [38]. In the *Drosophila* examples, imaging is predominantly at cell apices to capture the Zonula Adherens at which cortical actomyosin drives many cell behaviours, and at which E-Cadherin transmits tension.

The first step in putting numbers to morphogenesis is to quantify the rate of change of tissue shape at a spatial and temporal scale relevant to biological questions, described in sections (a) and (b) below. The details of the possible combination of cell behaviours that are responsible can then be quantified, as described in sections (c)–(f).

### (a) Strain rate basics

In one dimension (1D), the rate of change in length of a tissue is calculated across a time interval, typically between subsequent frames of a time-lapse movie. The strain rate, (where is the conventional symbol for a shape change or strain, and the dot for rate), is calculated as a change in length, scaled by the original length making it a unit-free proportional change, divided by the time interval, *δt*, to give a rate:
2.1where the log approximation can be used for small deformations (less than 20%). A positive (or negative) indicates a relative expansion (or contraction) rate. The length *l* can be calculated as the width of the tissue, or at smaller scales within the tissue. At the scale of distances between cell centres, for example, this will provide details on variation in the strain rate within a tissue, such as any biologically interesting gradient or pattern. In 1D, is a scalar and also a tensor, in that it captures the relative motion of points independent of any fixed reference frame (that is, the mean translation of points relative to the image coordinate system, or some landmark is not accounted for), and can be used to displace points and deform objects. Rearranging equation (2.1) gives
2.2

Thus the tensor changes the length of an object by the specified proportional rate. The same operation can be carried out to describe a change in 2D (or even three dimensions, 3D), but whereas is a single value in 1D, in 2D it becomes a 2 × 2 matrix (and in 3D, a 3 × 3 matrix), which can include not just length change but also rotation and shear movements.

For a 2D object or set of points that changes shape, could simply be phrased as a rate of change of area rather than length. However, how 2D tissue patches (domains) change shape is often not isotropic. Rather, they deform along a particular axis at a different rate or sign compared with how they change in the perpendicular orientation. That is, an elliptical description of change of shape, with independent and in perpendicular orientations, is more useful and complete. Higher-order descriptions of patch deformation are also possible, but it is not yet clear how relevant these will be to biological mechanisms. As well as changing shape, a 2D patch can also rotate. Hence in 2D, a complete linear description of the relative movements of a set of points in a plane requires four quantities (see figure 2 for example domains). These are (i) an orientation of the axis of greatest (principal) strain rate, (ii) the magnitude of that principal strain rate, , (iii) the magnitude of the strain rate perpendicular to the principal axis, and (iv) a rate of anticlockwise domain rotation. Strain rates are again calculated in proportional size change per unit of time, and rotation in radians per unit of time. Although these four quantities are sufficient to provide a complete linear description, there is a more elegant (albeit less intuitive) form in which these quantities can be precisely captured and handled mathematically, namely as a matrix tensor. What follows here is a brief explanation of how tensors work for describing tissue movements in 2D.

The two strain rates, orientation and rotation quantities are encoded in 2D in a 2 × 2 matrix tensor known as the tissue velocity gradient tensor, (pronounced ‘grad-v’). Note that each of the four elements of does not simply represent one of the four quantities above. Instead, the four quantities are mixed up across the elements of the matrix (see matrix in figure 2 examples), meaning that these quantities need to be extracted in specific ways (next section). However, has the elegant feature that average or cumulative shape change can simply be calculated by averaging or adding these matrices together, irrespective of any differences in orientation of the axes of deformation. Scaling up in space and time are therefore trivial. Strain rate data can, therefore, straightforwardly be presented versus time in one or a combination of ways—as spatially averaged instantaneous strain rates, as cumulative strains or as cumulative stretch ratios. can also be used to deform an object or displace a set of points in a similar fashion to equation (2.2). Note that an object that is not centred on [0,0] will also be translated by the transformation.

### (b) Calculating tissue strain rates

Some aspect of cells must be tracked, or registered, faithfully over time in order to measure the detailed change in shape of the tissue. Cell nuclei or cell centroids (the centres of mass, calculated from cell perimeter coordinates) are straightforward outputs of nuclear and cell membrane cell tracking, and both are ideal registration points. Even in the absence of cell registration points, for example, if cell tracking is problematic or as a precursor to cell tracking, alternative methods such as particle imaging velocimetry (PIV) and particle tracking velocimetry (PTV) can be used as a basis for calculating tissue strain rates [38].

All current strain rate methods start by calculating a very similar tissue velocity gradient tensor, before diverging in how they quantify and separate the contributions of component cell behaviours. A spatial and temporal scale at which to quantify strain rates must be chosen. The general approach should be to choose the minimal scale at which the component cell behaviours can be quantified and separated, since strain rates can easily be accumulated in time and averaged over space. The minimal scale at which cell rearrangement, cell loss or gain from the plane and cell division can be measured is four cells and their included cell–cell interfaces. Thus any domain size at or above this scale is applicable. Typically, the frame interval required to track cell shapes with automated methods is as short as that required to calculate strain rates that capture behavioural detail, since the speed of cell behaviours imposes a similar limit on both. Some methods require small size changes (less than approximately 20% strain between time points), hence care is needed in choosing methods and sampling interval. In principle, however, there is no upper limit to the sampling interval as long as cells can be faithfully matched between frames. Even in the absence of cell matching in time, carefully chosen multiple cell averages can be used to estimate cell behavioural strain rates [37]. The focus of most currently published methods has been on spatio-temporal domains of a cell and one corona of immediate neighbours upwards, and at time intervals such that the relative movement of cells is less than approximately one cell radius per frame.

There are three ways in which the tissue velocity gradient tensor, , is approximated in current methods. Note that it is not strictly that is calculated, but a discrete approximation to it, based on cell centroid movements rather than a continuous velocity field. As a result, nomenclature varies according to which approximation is used. In the first method, each of the four elements of the matrix **L _{T}** is calculated directly from a field of cell velocities. The gradients of the (linear) regression (i.e. best-fit line; figure 3

*a*) of velocity versus position are calculated in the following combinations [1,37,39]: 2.3where

**is the velocity in**

*u**x*and

**is the velocity in**

*v**y*. Similarly, PIV and PTV both generate velocity fields, the spatial derivatives of which are again the four elements of [38]. Note that the

*xy*locations and displacements of points in a domain are explicitly used to calculate velocity gradients.

The second approach to calculating is to map the evolution of ‘texture tensors' between time points. Texture tensors do not use *xy* locations explicitly, but instead describe the state of a domain by the distribution of lengths and orientations of links between the centres of all neighbouring cell pairs, irrespective of the locations of these links within the domain. Texture tensors quantify the average link length and the strength and orientation of any anisotropy in link lengths [40,41]. The evolution of the texture of persistent links (those not involved in a topological change) between time points is calculated as an optimized linear mapping, [35] (figure 3*b*, top line).

The minimal information needed to approximate in 2D are the coordinates of a deforming triangle. In the third so-called ‘triangle method’, therefore, triangles defined by the centroids of each trio of neighbouring cells are used. First, each triangle at each time point is completely described by a tensor that defines its shape and orientation relative to a universal reference triangle. Then, calculating a transformation tensor between the shape tensors in successive frames for a triangle, *m,* gives [29,42,43], which can be averaged over all triangles in a domain of cells (figure 3*b*, bottom line).

Information contained in is decoded in the following way (figure 3*c*). The rate of rotation is read from the off-diagonal elements of the ‘spin matrix’, Ω, calculated as
2.4where is the transpose of the matrix (i.e. with columns and rows swapped). The rotation-free strain rate tensor can be calculated as
2.5from which the magnitudes of the two perpendicular strain rates and their orientation can be extracted using standard functions in most numerical software packages. The strain rate magnitudes are often referred to in the language of matrices as ‘eigenvalues’, and their orientations as ‘eigenvectors’. Ω_{tissue} is known as the ‘antisymmetric’ part of , with the ‘symmetric’ part. Optionally, the strain rate tensor can be further separated into its trace (sum of the diagonal elements), representing the rate of area change, and a remaining pure shear or balanced convergence–extension motif (figure 3*c*). However, though this separation is commonly used, its use can be confusing. For example, if a stretch is imposed on a domain in one orientation by neighbouring tissues, a single positive rate of strain in that orientation would more closely reflect the biological mechanism than would an isotropic expansion added to a pure shear, though both are mathematically equivalent (figure 3*d*). The choice of representation should reflect the likely underlying mechanism.

Two aspects not accounted for in are the mean translation of the domain and any nonlinear or non-homogeneous behaviour within the domain. A mean translation velocity can be extracted from the regressions used to calculate elements of above. For example, the *x*-component of the translation velocity is the value of ** u** at the intercept where the best-fit regression line crosses the

*y*-axis in the plot in figure 3

*a*. See also the effect of adding translation to ‘cell’ trajectories in the example domains in figure 2. Non-homogeneous behaviour can be quantified in various ways, such as with higher-order (nonlinear) fits to the data, or from regression residuals (figure 3

*a*), or from the differences between cell coordinates at

*t*+

*δt*and

*t*after the latter have been deformed by . An example of cell behaviour that will produce zero tissue strain rates on average, but with strong unstructured residuals is tissue mixing, as can occur in endothelial sheets [44]. Structured residuals are also potentially biologically interesting, for example, if behaviour changes in a discrete rather than continuous manner across the tissue.

As approximated above, , the tissue velocity gradient tensor, describes the rate of deformation of a domain of tissue. In the next section, I will outline the various methods that are used to break this tensor down into the separate additive contributions of all the cell behaviours in the domain.

### (c) Behavioural strain rate decomposition

The goal of all strain rate decomposition methods is to quantify the separate additive contributions of all cell behaviours to the change in shape of the tissue: 2.6

Here is planar cell rearrangement, characterized by T1 processes, while is rearrangement into or out of the plane, through T2 processes, for example, due to the thinning of a multilayered epithelium or due to cell extrusion (which also captures apoptotic cell death since it is generally preceded by controlled extrusion [45]) (figure 4*a* and see Glossary of terms in box 1). I will make the assumption that domain rotation, as measured from cell centroids, Ω_{tissue}, is equal to the rotation of cells in the domain, Ω_{cell shape}, so that there is no rotation component to the other cell behaviours (but see Simple shear section).

The four published strain rate decomposition methods I will discuss can usefully be separated into two types. A ‘discrete method’, in two variants, makes direct use of topological changes and does not measure cell shape change directly. The discrete method has been applied to the pupal wing disc [29] and notum [35,40,47,48] and the chick primitive streak [38]. A ‘continuous method’ measures cell shape change directly, estimating the gradual process of intercalation by a difference method (see below). Topological changes are incidental to this method. It has been applied to the *Drosophila* germ-band [1,23,30], zebrafish trunk ectoderm [1] and the mammalian palatal epithelium [37] that elongate by a factor of 2 or more, and to the *Drosophila* amnioserosa [32,33] and salivary gland [34] tissues. Discrete methods are currently more comprehensive, while the continuous method currently lacks ways to separate out behaviour associated with each kind of topological change. However, I will sketch out below how progress might be made towards a comprehensive continuous decomposition. There are advantages to both kinds of method. Their suitability will depend on what tissues they are being applied to, what cell behaviours are present and what information is required.

For quantifying the contribution of particular cell behaviours to tissue deformation, information about cell shape and connectivity is required and this is normally extracted using (semi-)automated tracking of cell membrane fluorescence. The first behaviour to quantify is the cell shape strain rate, . Continuous methods measure cell shape change directly, finding the mapping that best fits the evolution of a cell's shape to the next image frame (figure 4*b*). This is performed on best-fit ellipse approximations to cell shapes (drawn as major and minor ellipse axes in figure 4*b*). Note that individual cell shapes can have an apparent rotation that is different from the domain tissue rotation, but all current methods assume that they do not (see Simple shear section).

Discrete methods do not measure cell shape strain rates directly. Instead, is calculated from all changes in the length and orientation of centroid links, or triangles of links, that are not involved in topological changes. That is, it contains both the effects of cell shape changes and the continuous sliding of cells past each other (intercalation) for a given frame interval, as long as these change the relative locations of cell centroids but do not involve connectivity changes. It is possible for cells to change shape and for this not to be detected by discrete strain rate methods. For example, if all cells in a domain elongate in one direction without changes in centroid location, area or connectivity, this behaviour will be cryptic (figure 4*c*). Thus the continuous and discrete methods capture different behaviours in . The next step is to quantify behaviour involving changes in connectivity.

### (d) Topological changes

Discrete methods excel at separating the remaining topological changes. They separately quantify three behaviours in consecutive image frames between which the topology changes. T1 processes are identified where cell connectivity flips between a tetrad of neighbours, T2 processes by the addition/loss of a new centroid and three associated cell–cell contacts, and cell division is identified by the approximate halving of a cell into two daughter cells, again with the addition of three new contacts (figure 4*a*). A deformation is attributed to each topological event in one of two ways (as for tissue strain rates, figure 3*b*). The first method calculates a strain rate tensor for each cell behaviour from the change in ‘texture’ of the subset of links connecting neighbouring cell centroids directly affected by each behaviour [35,38,40,47,48]. In the second method, the population of triangles connecting the cell centroids of trios of immediate neighbours that are directly involved in the topological changes are used [29]. Since T1 and T2 processes and cell division are generally rare events, discrete methods average over large domains (in space and/or time) in order to avoid producing noisy cell behaviour strain rates.

Continuous methods are used to calculate an intercalation strain rate tensor by subtracting the directly measured cell shape strain rate tensor from the tissue strain rate tensor (figure 4*d*) [1,37]. This captures the continuous process of cell rearrangement, to which topological changes are incidental. For example, during axis extension in the *Drosophila* germ-band, a decision seems to be taken some 5 min prior to neighbour exchange to start contracting the intervening cell–cell junction (between C and D in figure 4*e*), and relaxation of the length of the new junction after exchange takes a further 5 min. This continuous process of cell–cell slippage [1,46], in the middle of which cell topology changes, is captured by the continuous strain rate method. Indeed, comparing the continuous rate of intercalation rate during germ-band extension with the rate of topological change reveals a delay of some 5 min (figure 4*f*). Thus the continuous method correctly measures the onset of intercalation behaviour, which is significantly earlier than the onset of topological changes.

At the point of neighbour exchange involving two pairs of cells, all four cells meet at the centre of an ‘X’ cell–cell interface motif (at 10 min in figure 4*e*). In some tissues, neighbour exchange can involve more than four cells, which meet at a so-called rosette motif [49]. The cell rearrangements associated with rosette formation and resolution will be appropriately quantified as a local intercalation strain rate by each of the above kinematic methods. The local strength of the intercalation strain rate could be a useful proxy for rosettes in some situations or tissues, but the speed of rosette resolution is likely to be highly variable. Other analyses, such as identifying clusters of neighbour exchange events in space or identifying characteristic cell geometries through rosettes, could be used in addition to strain rate methods to highlight this particular behaviour.

Similar to T1 processes, T2 processes and cell division are more than just discrete topological events. Rather, they are ritualized behaviours that can take the focal cell and its surrounding cells some minutes to prepare for and to adjust to afterwards. Their influence on domain shape change could be calculated as some combination of cell shape change and cell rearrangement. That is, the continuous cell shape () and intercalation () strain rate tensors could be thought of as ‘atomic’ or ‘cassette-like’ continuous strain rate building blocks, out of which larger and more complicated behavioural classifications could be constructed. To do this, the expected start and stop times of these behaviours relative to topological changes would need to be defined, over which window net behaviour could be accumulated. Such an approach would provide a useful alternative to the discrete methods.

### (e) Continuous approach to cell division

I will now focus in more detail on how the net effect of cell division could be calculated. This is important because one might expect cell divisions that have a stereotypical orientation to drive tissue elongation. However, cell divisions oriented along the anterior–posterior axis of the zebrafish have been shown to be dispensable for axis extension [50] and body structures develop normally when cell division rates are altered [51]. Whether a mitotic event has a net contribution to the local planar shape of the tissue will depend on a number of factors, namely how the shape of the post-mitotic daughter cells combined differs from the pre-mitotic shape of the mother, how separated the daughter cells become, and how the local tissue accommodates the rounding up of the mother and then the separation of the daughters. This would therefore need to be assessed over the duration of mitosis, which from the start of prophase to cytokinesis is around 10 min in *Drosophila* epithelia [52], in addition to the post-cytokinesis time it takes for the daughter cells to ‘relax’ into the epithelium (figure 5*a*).

The discrete methods associate the tissue contribution of cell divisions instantaneously to the frame interval between which cell topology changes. The quantity that this yields depends on the precise geometry of local cell centroids. Hence a different frame interval might lead to a slightly different estimate, since the continuously evolving process of cell division will be imaged (or ‘frozen’) at different elapsed times from the division moment. The continuous methods currently quantify cell division either statistically from changes in cell number [37] or do not explicitly separate out its contribution [1], so here I will sketch out what the continuous approach to cell division would look like. The discrete methods could be extended to accumulate a net ‘cell division strain rate’ surrounding a cell division topological change, but the continuous approach seems better suited since it explicitly tracks exact cell shape changes. So how might this be done?

A net cell division strain rate will be extracted as an accumulation of some combination of the two atomic continuous strain rates, namely the cell shape and rearrangement strain rate tensors. Various biological parameters will need to be established to determine what to measure where and when. I will first consider the dividing cell only, ignoring its surrounding cells, and characterize strain rate contributions before and after the end of anaphase. This is the point at which the cytokinesis ring pinches a dumb-bell shape from the elongated anaphase cell (figure 5*a*). Up to the start of pinching, the cell shape changes of the mother cell will be accurately quantified in the cell shape strain rate tensor. From then on, there is a cell rearrangement component, which can be thought of as the two neighbouring cells adjacent to the cytokinesis ring invading into the cleavage plane [53]. This will be quantified by the intercalation strain rate tensor. Hence, if the moment when cells enter prophase can be identified, and cells can be tracked through cell division to the moment when daughter cells have equilibrated their cell shapes and locations relative to neighbours, an independent cell division strain rate tensor could be constructed as the sum of the cell shape strain rate tensors across the whole process, plus the intercalation strain rate tensors from the start of cytokinesis onwards.

However, how surrounding cells accommodate and respond to the above changes will ultimately determine whether there is any lasting contribution to tissue shape due to the mitosis event (ignoring subsequent interphase cell volume increase, or growth). If the local neighbouring cells temporarily and elastically accommodate the rounding of the mother cell and then locally flow to take the space between the separating daughter cells [54], then there will be no net contribution to the tissue strain rate tensor (figure 5*b*, left column). The cell mixing non-homogeneity will be measurable in the domain residuals. On the other hand, if the elongation of the mother cell area in the orientation of cell division results in similar movements of neighbouring cells, then the tissue domain will change shape as a result of the cell division (figure 5*b*, right column).

Thus a cell division strain rate tensor could indeed be constructed by integrating atomic continuous strain rates over a predefined time window from the start of prophase through to when daughter cells have relaxed into their new geometries, in the context of what the local neighbourhood does. Establishing the time window of cell division will need some understanding of the mechanisms driving cell division–associated movements, and of the compliance and fluidity of the local neighbourhood. A continuous approach to T2 processes (cells leaving or joining the epithelium) could be constructed along similar lines. In [37], an elegant alternative is proposed for when some 3D information is available.

### (f) Simple shear

The interpretation of domain rotation (extracted in equation (2.5)) is not straightforward. A combination of pure shear and rotation in a domain could be due to one, or a combination, of three mechanisms: pure shear (convergence–extension) and a global rotation of the domain, simple shear within the domain due to similarly oriented simple shear of cell shapes, or simple shear within the domain due to cells sliding relative to each other and not changing shape (figure 6). Quantifying and interpreting simple shear is a challenge for all decomposition methods due to the difficulty of extracting an independent measure of cell rotation. Continuous methods currently use the evolution of best-fit ellipses to cell shapes to calculate a cell shape strain rate. However, this shape information is incomplete in one aspect, namely that there is currently no way to identify registration points on cell membranes that can be tracked over time, with which to measure cell rotation [1]. Vertices where three cells meet are clearly useful points, but cell–cell interfaces appear and disappear as cells rearrange, divide, join or leave the epithelium. Given the fast rate of cytoskeletal turnover (half-life in the region of 30–60 s [55]), it is not clear whether subcellular long-term registration points in the cortex could or should be used. One could make the assumption that cell rotation is equal to the domain rotation, as both discrete and continuous [46] methods currently do, but there is no empirical evidence to justify this and this would rule out being able to distinguish cell shape from cell-sliding simple shears (figure 6*c,d*). A better understanding of the nature of vertices [56,57] and how cells slide past each other, or the use of a non-cortical marker of cell rotation might lead to practical solutions to this problem.

### Epithelial tilt and curvature.

Epithelia are rarely perfectly flat, being better described as quasi-2D surfaces [1] (figure 7). Hence, care is required if maximum intensity projections (MIPs) are used. The apices of some cell domains will be curved and tilted. Domain tilt is not an issue *per se* for strain rate calculations, which are calculated in proportional size change per unit time (a stretching domain has the same strain rate as a domain half the size stretching by half the amount), though estimates of velocity for tilted domains will be awry. If simply projected onto 2D, curved domains will be subject to nonlinear transformations, leading to errors in strain rate estimates. Errors can be reduced by flattening domains appropriately. Tilt and two orthogonal principal radii of curvature can be calculated across each domain from 3D points on the surface of the epithelium. The adjustment for a tilt of *θ* radians from the *xy* plane is to stretch the domain in the orientation of tilt by *x*′ = *x*/cos *θ*. For curved domains, these should then be stretched in each principal curvature orientation by *x*′ = sin^{−1}(*x*/*r*), where *r* is the radius of curvature, prior to strain rate analyses. Once calculated, accumulating or averaging strain rates across domains oriented in different planes remains an interesting unresolved problem. Other 3D issues remain to be fully incorporated into current methods. For example, the flow of an epithelium over a lip or corner will lead to transient shape changes as cells accommodate the acute curvature, for example by forming wedge shapes [63]. Similarly, a change in epithelial curvature will also introduce some kind of cell shape or arrangement change [64], which would need to be separated from other unrelated behaviours within the epithelium.

## 3. Discussion

Recent strain rate decomposition methods have begun to revolutionize how morphogenesis is quantified. Spatio-temporal strain rate maps in the various flavours of cell behaviour produced by these methods represent a significant advance in the precision and breadth of our description and quantification of phenotypes. Patterns revealed in these maps, such as spatially graded cell behaviours, their orientation and temporal sequence, suggest hypotheses about driving forces and constraints [29]. Mutant phenotypes might not be seen at the tissue level but become clear in the relative contributions of component behaviours. For example, during the first 30 min of *Drosophila* axis extension, the tissue is pulled from the posterior of the embryo by the posterior mid-gut invagination [30,58] at the same time as polarized Myosin II within the germ-band drives its convergence. AP-patterning mutants, in which the intrinsic convergence mechanism is impaired, still extend at wild-type rates. However, the relative contribution of cell rearrangement and cell stretch to tissue extension changes in favour of the latter, as passive cell stretch makes up for the lack of active cell intercalation [23].

Mechanical hypotheses based on strain rate patterns can then be tested by physical perturbation, such as laser ablation, by genetic manipulation and by optogenetic perturbation, and explored with mechanical modelling [29,30,37,48,59–61]. No one method on its own is likely to resolve a mechanical understanding, but complementary methods can reinforce each others' conclusions, allowing a mechanical understanding to be boot-strapped. There is also the tantalizing opportunity to develop a new generation of mechanical inference methods. Detailed strain rate maps can be combined with visual force sensing [62] or calibrated fluorescence of Myosin II motors as a proxy for active force generation [13] to infer tissue stress and mechanical properties.

Importantly, strain rate maps provide quantitative data on which to apply robust statistical comparisons. Much effort has been focused on developing ways to overlay multi-embryo data, coordinating embryos in space and time [23,29,32,33,35]. This leads to each *x, y,*(*z,*) *t* tissue location being associated with a population of strain rate values for each cell behaviour, from which within- and between-embryo variances can be calculated. Wild-type behaviour can then be compared to mutant behaviour, and phenotypic differences established in remarkable detail. However, the full statistical treatment of these multidimensional spatio-temporal maps remains to be established (and see box 2).

I have highlighted the differences between discrete and continuous approaches to quantifying the strain rate contributions of cell behaviours to tissue shape changes in domains of cells. Both use the relative movements of cell centroids, and then the discrete methods focus on changes in topology without directly measuring cell shapes, while continuous methods directly measure cell shapes to capture the continuous sliding of cells past each other. Discrete methods are elegant in their detail, but both methods have their advantages and drawbacks. The next stages in the development of continuous methods should focus on developing sets of rules to accumulate the net contributions of behaviours that involve topological changes. These will be built from the atomic continuous building blocks, namely cell shape and cell intercalation strain rates.

Domain shape changes emerge from subcellular mechanisms, notably cytoskeletal force generation, in the context of stresses and constraints from elsewhere [10]. In *Drosophila* embryonic epithelia, Myosin II contractility is organized into two subcellular populations, apico-medial and at junctions, and is often pulsatile [65–67]. This pulsatile contractility at the scale of single junctions and cell apices should have interesting signatures in the different flavours of strain rate. For example, apical contractions lead to junctional contraction and cell intercalation with interesting and quantifiable phase differences [58,68]. Understanding the subcellular mechanisms will inform how we should partition and interpret cell behaviour at the mesoscopic domain scale.

Unlike in many physical granular materials [69] where the grains are solid units, cells are active (often self-deforming), visco-elastically deformable units, with adhesions conferring a viscous friction on cell–cell contacts. Epithelia, nevertheless, have a cellular granularity that imposes constraints on tissue deformations. For example, the geometric arrangement of cells (grains) imposes constraints on what rearrangements are possible. It is impossible for cells to flow through each other (except in the case of cell division, which can facilitate this kind of fluid-like tissue behaviour) and simple shear needs a relatively straight boundary along which cells can slide past each other. The first step towards acknowledging geometric constraints would be to relate strain rates to the actual shapes and arrangements of cells that the strain rates are modifying. Cell shape change or area change can simply be related to the current shape and area of cells. More interestingly, it should be possible to classify the static arrangement states through which cells flow during cell intercalation, from which it should be possible to distinguish signatures of active cell-intrinsic from passive extrinsic mechanisms [70,71]. Thus far, no comprehensive theory links the static geometries and arrangements of epithelial cells to their kinematics, predicting constraints, but it now seems within reach.

## Competing interests

I have no competing interests.

## Funding

G.B.B. was supported by grant no. 15.23(k) from the Isaac Newton Trust and by Wellcome Trust grant no. 100329/Z/12/Z awarded to Prof. William Harris.

## Acknowledgements

I am grateful to Richard Adams, Alexandre Kabla, L. Mahadevan and Marko Popović for discussions on strain rate methods, and to Pedro Machado and Jocelyn Etienne for conversations about cell and tissue mechanics. I also thank the groups of Richard Adams, Bénédicte Sanson, Nicole Gorfinkiel, Alfonso Martinez-Arias and Katja Röper for collaboration on implementing strain rate methods, and thank Jeremy Green, Tara Finegan and Thomas Sharrock for comments on the manuscript.

## Footnotes

One contribution of 13 to a theme issue ‘Systems morphodynamics: understanding the development of tissue hardware’.

- Accepted December 5, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.