From heterogeneous morphogenetic fields to homogeneous regions as a step towards understanding complex tissue dynamics

During tissue development, morphogenesis is accompanied by cellular processes such as cell division, cell rearrangement, cell size and shape changes, apical constriction and apoptosis. The cellular processes are coordinated together, yielding collective cell migration and local deformation of each tissue region, resulting in convergent extension or epithelial folding. Furthermore, the local deformations of different tissue regions are coordinated too, resulting in large-scale tissue morphogenesis. Coordinations between invaginated mesoderm and covering ectoderm (Rauzi et al., 2015; Perez-Mockus et al., 2017) and between invaginated midgut and elongated germ-band (Collinet et al., 2015; Lye et al., 2015; Dicko et al., 2017) of the Drosophila embryo, between contracting wing hinge and expanding wing blade in Drosophila pupa (Etournay et al., 2015; Ray et al., 2015), or between invaginated neural plate and covering epidermal ectoderm of Xenopus embryo (Brodland et al., 2010), provide examples of how mechanical force generated in one region can drive large-scale deformation in adjacent regions. In these cases, the regions which behave differently are easily distinguished by specific gene expressions.

However, many tissues were found to be heterogeneous but without obvious boundary between such regions, leaving analysis limited to arbitrary regions drawn as a grid parallel to tissue axes, or regions expressing already known differentiation maker genes.

Measured tissue deformation rate showed a large heterogeneity (accompanied by a heterogeneity in cellular processes such as cell proliferation rate, cell division, cell rearrangement and change of cell shape) and smooth spatial variations across the tissue in Drosophila notum in a developing pupa (Bosveld et al., 2012; Guirao et al., 2015) (Fig. 1A,B), Drosophila wing blade (Etournay et al., 2015), blastopore lip of Xenopus gastrula (Feroze et al., 2015), chick gastrula (Rozbicki et al., 2015; Firmino et al., 2016), mouse palatal epithelia (Economou et al., 2013) and mouse limb bud ectoderm (Lau et al., 2015). Recent formalisms have enabled us to measure and link quantitatively cellular processes with tissue deformation (Blanchard et al., 2009; Guirao et al., 2015; Etournay et al., 2015; Merkel et al., 2017). These studies have shown that cellular quantities also vary smoothly across the tissue. In addition, the causal relationship between cellular processes and tissue deformation is not always trivial, making it difficult to identify regions that actively drive morphogenesis and that are passively deformed by adjacent regions.

To study the spatial regulation of morphogenesis at a tissue scale, we developed a new multi-technique pipeline to divide a tissue into sub-regions based on quantitative measurements of static or dynamic properties of cells or tissues. Our tissue segmentation pipeline consists of two steps and an optional third step: first, a fast tissue segmentation attempted several times with random seeding, then merging these multiple tissue segmentations into a single one, then, if necessary, smoothing the resulting region boundaries (Fig. 1C). In contrast to other image segmenting methods like the watershed algorithm, this method is designed to accept any kind of quantity of biological interest, not only a scalar but also a vector, a tensor, and combination thereof. Here, we have applied it to the morphogenesis of Drosophila pupa dorsal thorax and wing blade. They were divided based on the tissue deformation rate or how the cellular processes contribute to the tissue deformation. The sub-regions obtained showed distinctive patterns of deformation and cellular processes with higher homogeneity than those along tissue axes. Interestingly, the tissue segmentations based on the local tissue deformation rate and on the cellular processes included some similar regions, suggesting that the cellular processes were regulated similarly inside the regions, therefore resulting in homogeneous tissue deformations inside those regions.

Finding distinctive and homogeneous regions inside the heterogeneous tissue is achieved by segmenting the geometrical space while keeping the points inside each region as similar as possible to each other in the property space. Here, we call ‘property space’ any morphogenesis quantification measured in the tissue, whereas ‘geometrical space’ refers to the two-dimensional space of cell patch positions inside the tissue.

Given a set of objects, collecting similar objects to divide them into groups is generally a task of cluster analysis. However, the cell patches distribute both in the property space and geometrical space. On the assumption that expression patterns of genes responsible for morphogenesis make connected regions, and to study physical interactions between the regions, we aimed at getting connected regions. The initial tissue segmentation first defines a metric of similarity between cells, and then a tissue is divided into regions containing similar cells. A watershed algorithm is a widely used tool in general image segmentation and biology such as cell segmentation and segmentation of computed tomography data, but it requires the property space to be scalar and a clear boundary between regions. An alternative image segmentation tool, called region growing (Adams and Bischof, 1994; Ma et al., 2010) (Fig. 2A), was inspired by a study segmenting mouse heart based on cell polarity (Le Garrec et al., 2013).

To validate the algorithm, we first tested segmentation on a simple example, namely the change in cell patch areas from 12 to 32 h after puparium formation (APF) (Fig. 2B). The overall change in cell patch areas defines the total tissue growth, whereas spatially heterogeneous changes in cell patch areas result in local deformation, changes in tissue region proportions and overall tissue shape change. Technically speaking, the change in cell patch areas is a scalar field, defined as the trace of the tissue deformation rate tensor. The region growing succeeded in finding expanding regions in posterior, lateral posterior and lateral parts, and a shrinking region in the anterior part.

However, the results varied dependent on the initial seeds. In contrast to a segmentation of immunostained images, where a true segmentation is well defined, the morphogenetic properties vary continuously with space, making it difficult to determine and validate the resultant segmentations. The silhouette, a measurement of region homogeneity (the silhouette of an object would be 1 if it was similar to all objects in the same cluster, and −1 if it was more similar to objects in other clusters), differed from one segmentation to the other (Fig. 2C). To assess the significance of the homogeneity, we compared it with the average silhouette of randomly made control segmentations. Some of the region growing results had a low silhouette, even lower than that of half of the control segmentations (Fig. 2C), which means they were lacking any significance.

Because of the random initial seeding, we do not know which results should be compared with gene expression patterns or fed forward to a study of mechanical interactions between the regions. For practical applications, we need a single segmentation result for a given morphogenetic property.

To obtain a single tissue segmentation, we turned to consensus clusterings. In fact, as resultant segmentations of the region growing were dependent on randomly given initial values, we ran multiple trials and merged multiple segmentation results into a single one. Given multiple partitions, the consensus clustering returns the partition which is the most similar to all of the initial partitions. We tried several consensus clustering algorithms, and found the ‘label propagation on a consensus matrix’ (Lancichinetti and Fortunato, 2012; Raghavan et al., 2007) returning regions similar to the results of the region growing.

The label propagation on a consensus matrix converted multiple tissue segmentations into a weighted graph in which weight of an edge represented a frequency of segmentations in which incident vertices (points) belonged to the same region (Fig. 3A). Then labels on the vertices were propagated according to the weight so that the same label was assigned to points which were frequently included in the same region among the given multiple region growing segmentations.

The label propagation returned results similar to the region growing segmentations (Figs 2B and 3B). Also, the label propagation results were more similar to each other than results of region growing, assessed with adjusted Rand indices (ARI), a measurement of similarity between two partitions (ARI of identical partitions would be 1). ARI were 0.50±0.21 (mean±s.d.) among the results of the region growing and 0.97±0.02 among the results of the label propagation. They showed similar average silhouette values, similar to median of those of region growing results, but smaller than the highest value of the region growing results (Fig. 3C). The average silhouette of the label propagation result was higher than those of 99.95% of the randomly made control segmentations.

However, a consensus clustering algorithm ignores original properties of objects in principle and divides the objects only based on how they were divided among given partitions, and thus it might return disconnected regions and a zigzag boundary between them. Some segmentations in Fig. 3B also included disconnected regions, marked in gray.

For the cases of a complex boundary and disconnected points, we prepared an optional step to smooth the boundary and remove disconnected points when needed. To smooth the consensus region boundaries, we employed the cellular Potts model, which simulates dynamics of a cellular tissue by calculating energies of cells from their geometry, trying to decrease the total energy. In our application for boundary smoothing, the energy was lower when the region boundary was shorter and the homogeneity was higher (Fig. 4A). Resultant regions were evaluated by homogeneity and circularity (Bosveld et al., 2016), which represents the smoothness of the boundary. Parameters for the cellular Potts model were screened so that the circularity was higher than a given value and the homogeneity was as high as possible. The result was not affected greatly by extending the duration of the simulation, and thus we stopped the simulation when it was appropriately smoothed.

This smoothed boundaries and removed disconnected cell patches (Fig. 4B,C) while keeping the average silhouette value higher than those of 99.5% of the randomly made control segmentations (Fig. 4D). As the cellular Potts model implementation includes the Metropolis update, i.e. choosing a pixel randomly and updating the pixel by probability according to a change of the energy, resultant smoothed segmentations varied among different trials even with the same parameters and initial segmentation. We therefore iterated the cellular Potts model smoothing 50 times and integrated its results using the label propagation algorithm again.

This gave a pipeline of the region growing, the label propagation and the optional cellular Potts model to divide a field of property (scalar, tensor, or any kind of value with metric) into regions. The resultant regions are homogeneous, where points in each region are more similar to each other than to points in other regions.

We tested whether our method can segment a simulated tissue using a normal cellular Potts model. We prepared a 2D tissue model in which two types of cells with different surface tensions were assembled (Fig. 5A). The tissue was geometrically compressed in the horizontal direction and extended in the vertical direction, so that the cells retained their area but they are vertically elongated (Fig. 5B). Because of the surface tension, the cells tried to minimize their perimeter, and thus they were rearranged and rounded (Fig. 5C; Movie 1). The tissue was split in a 4×7 grid, and the deformation and cellular processes after the compression were measured in each cell patch. The measured deformation and cellular processes were averaged among 12 simulations. The cells were rearranged similarly between the two types of cells (Fig. 5D; Movie 1). However, when we tried the segmentation based on the cell rearrangement, the cells were successfully distinguished (Fig. 5E). When the tissue was segmented into three regions, it included a small and disconnected region after the first label propagation and smoothing (Fig. 5F,G). The cell rearrangement rate was slightly higher among the cells with lower surface tension (Fig. 5H). In general, lower surface tension would allow larger fluctuation, and the larger fluctuation might have facilitated the faster cell rearrangement. The two regions obtained showed a significantly high homogeneity (Fig. 5I). With these results, we confirmed that our method could properly identify groups of cells with different mechanical properties based on their behavior.

We now turn to property spaces better representing tissue morphogenesis. In Guirao et al. (2015), tissue deformation rate (G) and underlying cellular processes, cell division (D), cell rearrangement (R), cell shape change (S) and cell delamination (A) were quantified into tensors. The tensors were obtained from the change of texture averaged over 20 h from 12 h APF to 32 h APF or over 2 h at each time point. By comparing the tensors, for example, one can check whether cell divisions and cell rearrangements elongated tissue in the same direction or attenuated each other. In the same way, by comparing the tensors of deformation rates with a unit tensor which has the same direction of elongation as tissue deformation rate, one can estimate an amplitude of the tissue deformation rate and how much the cellular processes contribute to the tissue deformation in both terms of contraction/expansion (isotropic deformation) and narrowing/elongation (anisotropic deformation) (Guirao et al., 2015). They are scalar values and denoted by G_// for the tissue deformation rate, D_// for cell division, R_// for cell rearrangement, S_// for cell shape change and A_// for cell delamination. For the sake of clarity, we call the tissue deformation rate and the cellular processes averaged over the whole 20 h from 12 to 32 h APF ‘time-average’ tissue deformation rate and cellular processes.

The effective contributions averaged over the whole tissue showed dynamic time evolution (Fig. S1), with a large peak of cell division and cell shape change at around 16 h APF, a second small wave of cell division at around 22 h APF, and a gradual increase of cell shape change and cell rearrangement. The effective contributions also showed large variance across the tissue at each time point. Therefore, we included the time evolution in the property space. Assuming that there are two regions in a tissue where the tissue expands, the first region expands during 14-17 h APF and the second region expands during 25-28 h APF, resulting in similar size changes, then the two regions cannot be distinguished by the time-average expansion rate. To distinguish them, we compared a property at each time point and summed up its difference through the whole time. When two cell patches always behaved similarly, then the difference at each time point is small and so the total difference is small too, whereas cell patches with deformations occurring at different times are separated at each time point and thus the total difference gets large. In contrast with time-average, we call the sum of difference at each time point ‘time-evolution’.

We first divided the tissue based on time-average and time-evolution of tissue deformation rate. The similarity was given by Euclidean distance of tensors. The notum was divided into anterior-middle-posterior and medial-lateral regions by both the time-average and time-evolution, whereas the middle regions were smaller and the middle lateral region extended medially in the segmentation based on time-evolution (Fig. 6A-D).

Next, we divided the tissue based on time-average and time-evolution of cellular processes. The amplitude of tissue deformation rate and effective contributions of cellular processes were combined in a vector, and their similarity was given by Euclidean distance of vector. In contrast to the segmentations based on the time-average and time-evolution of tissue deformation rate, the segmentations based on time-average and time-evolution of the cellular processes were dissimilar to each other (Fig. 6E-H). The segmentation based on time-evolution of cellular processes included a posterior region, a large anterior region, a neck-notum boundary region, lateral posterior region, a middle boundary region and a lateral region (Fig. 6H).

Although a change in the minimum circularity for the boundary smoothing did not affect the segmentation based on the time-average tissue deformation rate (Fig. S2A), the shape of the boundary changed dependent on the minimum circularity for the segmentations based on the other three property spaces (Fig. S2B-D). The cellular Potts model returns regions with circularity higher than the given minimum value if it is the most homogeneous segmentation.

We also tried dividing the tissue into various numbers of regions (Fig. S3). In many cases, an increase of the number of the regions resulted in subdividing the regions already obtained. In addition, when the number was too large, some results of the first label propagation included small regions similar to one observed in the segmentation of simulated data (Fig. 5F). Those small regions were absorbed into surrounding regions during the smoothing by cellular Potts model. Then the final label propagation tried to integrate regions smaller than the final segmentation, returning small and sometimes disconnected regions (Fig. S3, third column, from third to bottom rows; fourth column, sixth and bottom rows), again similar to the simulated data (Fig. 5G). Thus, the existence of small regions suggest that it is oversegmented.

Both of the segmentations based on time-evolution of tissue deformation rate and effective contributions of cellular processes included the large anterior region, the middle boundary region, the lateral posterior region and the posterior region, although the anterior and posterior regions were divided into medial and lateral sub-regions in the segmentation based on the tissue deformation rate. Fig. 7A-D shows overlap between segmentation based on time-evolution of cellular processes (Fig. 6H) and the segmentations based on time-average tissue deformation rate, time-average cellular processes and time-evolution of cellular processes (Fig. 6B,D,F) or a conventional large grid parallel to tissue axes. The middle lateral and posterior lateral regions in the segmentation based on time-evolution of tissue deformation rate, and the middle boundary region and lateral posterior regions in the segmentation based on time-evolution of cellular processes, overlapped each other (Fig. 7C). We also evaluated the overlap between the segmentations by ARI (Fig. 7E). Despite the difference between the anterior sub-regions, the segmentations based on time-evolution of tissue deformation rate and cellular processes overlapped each other more than the others.

Next, we evaluated the homogeneity of the obtained regions. The time-evolution of tissue deformation rate was similar among cells inside regions of the segmentations based on time-average and time-evolution of tissue deformation rate, except the middle-lateral region identified based on the time-average tissue deformation rate (Fig. 8A,B). On the other hand, the large grid segmentation showed large heterogeneity in the posterior regions (Fig. 8C). The average silhouette value of the segmentation based on the time-evolution of deformation rate was higher than that of 99.5% of the control segmentations (mean silhouette for label propagation, 0.0568; for smoothed regions, 0.0600; the maximum of the smallest 99.5% of control mean silhouettes, 0.0425) (Fig. 8D). The mean silhouette of the segmentation based on time-averaged tissue deformation rate was also higher than 95% of the control segmentations (mean silhouette for label propagation, 0.0336; for smoothed regions, 0.0260; the maximum of the smallest 95% of control mean silhouettes, −0.0098). However, the mean silhouette of the conventional grid segmentation was close to the median of the control segmentations (mean silhouette, −0.0694).

Also, the time-evolution of cellular processes was homogeneous inside the regions of the segmentation based on time-evolution of cellular processes, but not in segmentation based on time-average of cellular processes nor in the grid (Fig. 8E-G). The mean silhouette value of segmentation based on time-evolution was higher than 99.995% of control segmentations (mean silhouette for label propagation, 0.147; for smoothed regions, 0.1448; the maximum of the smallest 99.995% of control mean silhouettes, 0.124), whereas that of segmentation based on time-average was smaller than 5% of control segmentations (mean silhouette for label propagation, 0.0125; for smoothed regions, −0.0385; the maximum of the smallest 95% of control mean silhouettes: 0.0241) (Fig. 8H).

Our tissue segmentation is designed to divide a tissue into regions homogeneous in a given property space, and the homogeneity of either tissue deformation rate or cellular processes in the segmentations based on each property demonstrated that the pipeline worked (Fig. 8B,F). However, it does not ensure the homogeneity of the regions in other property spaces that are different from the property space on which our segmentation was performed. Fig. S4 shows heat maps of silhouette values measured in different property spaces. Even though the homogeneity in the regions differed among the different property spaces, the segmentations based on time-evolution of tissue deformation rate and cellular processes showed higher homogeneity than the others also in the property spaces of deformation rates owing to cell divisions, cell rearrangements and cell shape changes.

We projected the regions divided based on the time-evolution of cellular processes onto the actual cell map, and found that the anterior and posterior regions corresponded to scutum and scutellum, and the middle boundary and lateral posterior regions covered the scutum-scutellum boundary (Fig. 9). This result demonstrates that the obtained regions corresponded to the anatomical features, and cells were behaving differently between the anatomical regions. Note that the segmentations based on time-averaged tissue deformation rate or cellular processes did not match the anatomical features, indicating that the cells in the anatomical regions are actually regulated temporally.

Fig. 10 shows plots of the average effective contributions of cellular processes in each region of the segmentation based on time-evolution of cellular processes. The second peak of cell division was observed only in the posterior regions and the middle boundary region, consistent with previous studies providing maps of the number of cell divisions (Bosveld et al., 2012; Guirao et al., 2015); however, we also found a small first peak of cell division in the lateral posterior region. Plots of average cellular processes also differed from each other among regions in the segmentation based on tissue deformation rate, but this was less distinctive than in the large grid segmentation (Fig. S5). Distances between the plots in Fig. 10 were 0.65±0.16 and those for the segmentation based on tissue deformation rate were 0.63±0.20, larger than those for the large grid segmentation (0.44±0.14). This result demonstrates that cellular processes in the obtained segmentations were more distinctive than those in the conventional grid.

To demonstrate the generality of our method in dividing a tissue, we performed the same segmentation and analysis in the Drosophila pupa wing blade. During 15-32 h APF, the wing blade is elongated in the proximal-distal direction by a contracting wing hinge connected to the wing blade proximal side, whereas its distal side is anchored to the cuticle via Dumpy (Etournay et al., 2015; Ray et al., 2015). The wing hinge contraction also narrows it in the anterior-posterior direction and induces shear strain in wing blade proximal anterior and posterior regions (Fig. 11A,B). We performed tissue segmentation for the wing blade based on time-evolution of tissue deformation rate (Fig. 11C) and cellular processes (Fig. 11E), which divided it into four regions. In both cases, the wing blade was divided into anterior, middle, posterior and distal regions. All regions showed positive silhouette values (Fig. 11D,F), and their averages were significantly higher than the average silhouette values of control segmentations (Fig. 11G,H). Like the notum, dividing the wing blade into a larger number of regions also subdivided already obtained regions (Fig. S3). Plots of the effective contributions of cellular processes also showed distinctive patterns between the regions: the cell division showed a small contribution in the anterior region, the cell rearrangements dominated the tissue deformation around 26 h APF in the anterior and posterior regions, and the cell shape changes showed two peaks around 17 and 22 h APF in the distal region (Fig. 11I). Projection of the four regions onto the cells showed a difference between the regions and interveins – the posterior region roughly corresponded to the proximal posterior intervein and the boundary between the anterior and distal regions corresponded to the L3 vein (Fig. S6).

This study demonstrates that the pipeline of the region growing, the label propagation on the consensus matrix and the boundary smoothing by cellular Potts model could divide a deforming heterogeneous tissue into homogeneous regions based on any prescribed quantity. Using this segmentation method, we divided the developing dorsal thorax and wing of Drosophila pupa based on their morphogenesis, and found regions with a distinctive tissue deformation rate and underlying cellular processes.

The tissue segmentation based on morphogenesis differs from conventional image segmentation and cell segmentation algorithms in terms of quantity and supervision. Marking a human in a picture, an organ in a histological image or dividing cells in a microscopic image had been done manually, and thus results of the watershed algorithm or the artificial neural network could be supervised and corrected with the manual segmentations. On the other hand, the segmentation based on morphogenesis is hardly ever done manually for several reasons. First, the morphogenesis was quantified as multiple tensor fields with time evolution, and thus it is hard to visualize them in a 2D image for manual segmentation. Second, it is not easy to evaluate whether a given region actually corresponds to genetical/mechanical regulation of morphogenesis. Therefore we looked for a method which divides a tissue based on any prescribed quantity and returns regions with smooth boundaries. Region growing is a conventional and simple method of image segmentation, and only requires a property space to be metric. The varying results of the region growing were given to the label propagation and cellular Potts model to produce a single tissue segmentation, and the result was evaluated by region homogeneity boundary smoothness.

The notum segmentations based on time-evolution of tissue deformation rate and effective contributions of cellular processes returned similar regions corresponding to the scutum, scutellum and the boundary between them. Using the tissue deformation rate, the scutum and scutellum regions were divided into medial and lateral sub-regions. As the vector of effective contributions ignores the direction of deformation, the two sub-regions could be interpreted as regions of similar underlying cellular processes but deforming in different directions. However, the middle boundary region and the lateral region given by the cellular processes both overlapped with the middle boundary region given by the tissue deformation rate, and could be interpreted as regions with different cellular processes but of similar tissue deformations. The wing blade was divided into anterior, middle, posterior and distal regions based on both the tissue deformation rate and cellular processes, but the regions did not match the wing veins pattern.

Silhouette analysis showed that the segmentations based on time-evolution of deformation rate and cellular processes included regions homogeneous in various property spaces, whereas the conventional grid segmentation included heterogeneous regions. This method has some limitations. It cannot segment a small region; but such region might disappear during the boundary smoothing, resulting in disconnected regions in a final segmentation. Also, it is hard to determine the number of regions. We only know that a tissue might be oversegmented when it includes the small region. In a practical application, the tissue shall be segmented into various numbers of regions with various minimum circularities, and one of them can be chosen by comparing it to gene expression patterns, analyzing cell behaviors inside each region or other characteristics of interest.

In conclusion, we built a method to divide a tissue based on any prescribed property space. This allows the study of spatial regulation of various processes, where the property space should be chosen for the process of interest. For example, to study spatial regulation of cell division orientation, the property space may be prepared from, instead of the local deformation rate and cellular processes, the tensor field of cell division and known regulating factors such as cell shape, localization of planar cell polarity proteins and tension on cell-cell interface, and then resultant regions can be compared with gene expression patterns. Also, this method is not dependent on how the morphogenesis was quantified, and one can include rotational movement by anti-symmetric strain rate tensors, or 3D deformation by using voxels instead of pixels.

This section describes in detail the quantification tools. Morphogenesis data result from the quantification of local tissue deformation rate and underlying cellular processes as described in Guirao et al. (2015). The similarity between two tensors is quantified by the standard Euclidean metric. The homogeneity of a quantity within a given region, i.e. the similarity between measurements of this quantity within a region, is measured by silhouette, a standard tool of cluster analysis. For a measurement of similarity between tissue segmentations, we use the Rand index, which indicates how well two data clusterings agree.

Quantification of local tissue deformation and underlying cellular processes was performed as previously described (Guirao et al., 2015). Briefly, Drosophila nota expressing GFP-tagged E-cadherin were imaged. The notum movies were split in a grid (with patches about 20 μm wide) at the first frame (Fig. 1C), 12 h after pupa formation (APF). The local deformation rate and the cellular processes were measured in each cell patch through the development, as described below.

Epithelial cell contours were detected automatically using a watershed algorithm, cells were tracked, adjacencies between cells were listed and relative positions of adjacent cell centers were recorded. The tissue deformation rate, denoted by the symmetric tensor G, was obtained from changes of relative positions between neighbor cells over 20 h from 12 h APF to 32 h APF, or over 2 h at each time point when recording the time evolution. The tissue deformation rate G was then decomposed into cell shape change S and deformation accompanied by change of cell adjacency, which was further divided into cell division D, cell rearrangement R and cell delamination A, which are symmetric tensors too.

Similarity of morphogenesis between different cell patches was defined as follows. For expansion/contraction of an area (isotropic deformation), similarity was given by difference in expansion/contraction rates.

We took the average silhouette value over all points (cell patches) as a measurement of homogeneity of a given segmentation. For a significance test, tissue was segmented randomly (see below) 20,000 times into a given number, and we got thresholds above which the highest 5%, 0.5% or 0.05% of the average silhouettes were found. The average silhouette of given regions was compared with those of the control segmentations with the same number of regions.

For a measurement of similarity between tissue segmentations, we use the Rand index, which indicates how well two data clusterings agree; its value is 0 if the clusterings entirely disagree and 1 if they entirely agree. Its corrected-for-chance version, the ARI, is a more meaningful quantity: it is the Rand index compared with its value expected for the random case, and its value can be negative.

The pipeline was implemented by custom Matlab scripts, in three steps (Fig. 1C). The Matlab scripts are available at GitHub (https://doi.org/10.5281/zenodo.4270726). Also, see supplementary Materials and Methods for pseudo codes of the scripts.

The initial tissue segmentation first defines a metric of similarity between cells, and then a tissue is divided into regions containing similar cells. This approach was inspired by a study segmenting mouse heart based on cell polarity (Le Garrec et al., 2013). On the assumption that expression patterns of genes responsible for morphogenesis make connected regions, and to study physical interactions between the regions, we aimed at getting connected regions.

The algorithm ‘Region growing’ (Adams and Bischof, 1994; Ma et al., 2010) is an image segmentation method using a process similar to k-means clustering, starting from randomly given seeds (corresponding to ‘means’ in k-means clustering), segmenting an image with the seeds followed by update of the seeds within the regions, and iterating this process until convergence (Fig. 2A). The tissue segmentation is carried out by growing regions from the seeds collecting pixels adjacent to the growing regions, and so the resultant regions are connected.

Initial seeds were randomly chosen from data, and regions were expanded by adding a pixel (cell patch) adjacent to a region and the most similar to the seed of the region in the property space one by one until all pixels were assigned to one of the regions. The seeds were updated to pixels closest to centroids of the regions, averages of the regions in the property space were given as a property of the seeds, and then regions were expanded again from the seeds. These region expansions and seed updates were iterated until convergence was reached.

The control segmentations were made using an algorithm similar to the region growing but ignoring the similarity between points. From randomly given seeds, regions were expanded by adding a pixel adjacent to a region, where the added pixel was chosen randomly from all adjacent pixels, until all pixels were assigned to one of the regions. Therefore any obtained region is connected.

To merge multiple segmentation results into a single one independent on the metric, we used a label propagation algorithm on a consensus matrix, which takes multiple partitions and returns a consensus partition which is the most similar to all partitions (Lancichinetti and Fortunato, 2012; Raghavan et al., 2007).

For a division of n points, 50 independent trials of region growing were converted to a consensus matrix, the entry of which at i-th row and j-th column indicates a frequency of partitions in which i-th point and j-th point were in the same cluster. The entries lower than a given threshold were set to 0. The label propagation started by assigning a different label to each point. Then the label of randomly chosen i-th point was updated to one that was the most weighted by the consensus matrix, where ij element gave the weight to a label of j-th point. The label update was iterated until convergence. The threshold for the consensus matrix was scanned between 20 and 80% so that a resultant partition contained the same number of regions as the initial partitions.

When updating the label for a randomly selected pixel a, a target label was randomly selected from neighbors of a, and then Hamiltonian change was calculated. The label of a was updated to the target label with probability ⁠, where denotes change of by the change of label of a, and T is the fluctuation allowance. In the present application, the updates of labels were iterated 50 times. For resultant regions, a circularity C was calculated, where it was defined as C=4π×area/perimeter² (Bosveld et al., 2016). The parameters J, λ, h and T were screened for resultant regions with the highest homogeneity and circularity larger than a given threshold value. The screening was performed in a pairwise testing manner on a grid, and the grid was converged to the highest homogeneity. With the screened parameters, the boundary smoothing was iterated 50 times, and the results were integrated again by the label propagation on a consensus matrix algorithm.

For the simulation of cells in a compressed tissue, the cellular Potts model was implemented by custom Matlab scripts, which are available at GitHub (http://doi.org/10.5281/zenodo.5016684). It is simulated on a torus surface so the cells made a periodic pattern. The initial configuration was prepared by the Voronoi tessellation from 600 randomly scattered points. The points were randomly selected from 864×150 pt² plane in a sequential manner so that all points were at least 10 pt away from the others. Cells on the left side were assigned the low line energy J_soft=1 and cells on the right side were assigned the high line energy J_stiff=4. All cells were assigned the same compressibility λ=1. For the compression, the initial configuration was transformed into a 480×270 image. Then dynamics of the cells was simulated for 975,000 updates.

To be processed by the tools developed in Guirao et al. (2015), images of the cells were projected on a plane and converted to cell segmentation. As all pixels belonged to one of the cells in the cellular Potts model and the cells were not separated by boundary (watershed) pixels, we first labeled the pixels with ‘perimeter’ or ‘inner’, where the perimeter pixels were adjacent to different cells and the inner pixels were enclosed by the perimeter pixels. Then the perimeter pixels were re-labeled one-by-one to inner so that all inner pixels were connected in each cell, all perimeter pixels were adjacent to inner pixels of the same cell, and no inner pixel was adjacent to inner pixel in different cells. The labels were updated as much as possible, and remaining perimeter pixels were taken as the boundary.

The authors thank Floris Bosveld for imaging, Stéphane U. Rigaud for cell segmentation and tracking, Yohanns Bellaïche for useful advice, and all team members for insightful discussions.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.199034.