Image analysis of self-organized multicellular patterns

1 Introduction

Spatial organization of tissue and multicellular pattern formation are crucial biological processes in development, wound repair, and maintenance of tissue [1]. Understanding the underlying mechanism is of general interest and may be helpful for generation of organoid cultures or in tissue engineering. Tissue organizational processes are considered to be regulated by various physical and chemical factors, such as differently diffusive morphogens, cell–cell contacts, microenvironments and cell mechanics [2]. Accordingly, various theoretical models, often based on reaction-diffusion models, try to explain underlying principles.

In most tissues cells are emended in a complex extracellular matrix, a mesh of fibrous polymers attached to each other and to cells by a series of complex molecular links. The cells organize the geometrical and physical topography as well as the mechanical properties of the matrix, and are in turn responsive to these properties. It is commonly believed that the sensitivity of cells to the mechanical properties of their environment depends on their ability to exert contractile forces into their surroundings [3]. Indeed, many cell types including fibroblasts generate contractile forces via the actin cytoskeleton. The resulting matrix deformations generated by the cells may be felt by other cells at long distances and thus may lead to elastic interactions of cells [4], [5]. An interesting consequence of these interactions can be a spontaneous alignment and pattering of cells. Despite the implications of this organizing principle the pattering phenomenon has not been extensively studied in cell mechanics experiments. Interestingly, patterning phenomena can also be observed in simple 2D in vitro cultures (Figure 1), a fact well known by cell biologists. In order to investigate systematically how mechanical interactions of cells may govern their self-organization in a soft elastic environment we developed compliant elastomer substrates for in-vitro studies. On such functionalized surfaces we studied temporal aspects of pattern formation observed by long-term microscopy (phase contrast). However, the analysis of cell position and alignment direction requires either manual analysis or the application of specifically automated image analysis.

Figure 1

A culture of dermal fibroblasts on a collagen-coated compliant (70 kPa) PDMS-substrate at two different time points. On the left image (day 2 of culture) cells are at a low density and no apparent order is visible. On the right image (day 7 of culture), cells are in full confluence arrange in a pattern with high orientational order.

In our approach the overall pattern is analyzed based on the individual orientation of single cells. The required explicit identification and measurement of biomedical objects must respect heterogeneous representations of characteristic features such as shape or texture of a single type of object. In addition, automated processing of large series of images must also cope with technical limitations and artifacts such as inhomogeneous densities of objects or illumination within single exposures as well as in the series of images. In this approach a three-step multi-scale object oriented approach is proposed to extract visually plausible regions, identify relevant objects and measure the orientation and center of mass. This ensures a comprehensible, reproducible and therefore quantitatively assessable automation.

2 Cell experiments

3 Automated multi-scale object analysis

Cells are identified in three steps in digitized micrographs (Figure 2) [7]. At first an image is partitioned on different scales by a causal morphological scale space analysis, where the edges represent the zeroes of the second derivatives (Figure 3). This corresponds to human visual perception and is independent of local variations in illumination. In the second step the hierarchical structure of the regions is transferred into a tree where each node holds a set of attributes describing the region formally such as mean grey value, size, curvature, texture, etc. (Figure 4). In the third step regions that represent cells are extracted from the tree by application of a task specific top-down query on the attributes of the tree nodes.

Figure 2

A multi-scale representation of how an image (left) is transferred into an attributed tree representing the region hierarchy (middle). Relevant regions are extracted by specific attribute queries (right).

Figure 3

Multiscale partitioning of a phase contrast micrograph of 800 × 600 pixels (top left) with scale 1 (top right) scale 3 (bottom left) scale 5 (bottom right), where scale corresponds to a spherical template with radius 1, 3 and 5 pixels, respectively.

The task specific query is generated by heuristic analysis of few images and then applied to the entire series (Figure 5, Rules 1, and 2). On the extracted regions the major axis angle and the centre of mass describe the regions position and orientation (Figure 6). By considering the ratio of minor/major axis the general elongation is characterized. The process is reproducibly applicable to arbitrary series of images.

By inspection of some regions representing cells in the sample image (Figure 4) the following simple query has been composed in the domain specific language (DSL) of the multiscale analysis software (www.morphoscope.com).

Figure 4

In the hierarchy browser the attributes of single regions representing cells (left) can be inspected manually (right).

Rule 1

DECLS:

conv := (4 * area * 3.14176)/(pow((perimeter),2.0));

#DECLS

RULES:

TD:

(meangrey < 90) && (area < 5500) && (area > 500) && (conv < 0.2)

#RULES

In this case meangrey denotes the mean grey value of the corresponding region between 0 and 255, area is the absolute number of pixels and conv represents the roundness which is computed on-the-fly by the ratio of the regions squared outline in pixels (perimeter) and the area proportionally to pi. The declaration of conv is also given in the DSL (DECLS: Section). The extraction rule is defined in the RULES: section and applied in a top-down scheme (TD:).

Since the rule is based on heuristic and arbitrary analysis of single regions it is not a single and sound solution. For instance the following rule yields the same result as Rule 1 on the sample image (Figure 5).

Figure 5

Application of Rule 1 and Rule 2 yield the same result when applied to the sample image.

Rule 2

RULES:

TD:

(meangrey < 90) (area > 1000) && (greyentropy > 0.80)

#RULES

Here the attribute greyentropy denotes the entropy of the grey value histogram of the region. The similar results of two different rules indicate a certain degree of stability of the heuristic approach yet they underline the arbitrariness of the rule based object extraction. However, the blobs of identified regions show a visually comprehensible representation of the orientation distribution of the single cells (Figure 6).

Besides the five attributes for representation of the blobs any set of region features may be defined in the DSL and provided for further analysis such as orientation distribution.

This approach of region extraction by separation of image analysis and content description allows for rapid verification of different extraction queries without the need of parameterization of low-level image processing techniques, or low level development of new algorithms. The knowledge of the biomedical expert is represented by the query rules.

In the image (Figures 3–6) 96 regions were identified by Rule 1. 5 Regions were background (false positive) in 3 cases a single cell was split into two regions. However, a systematic validation of the results found by the algorithm against a manually generated ground-truth is needed, to estimate error rates of automated high-throughput screenings.

Figure 6

Visual representation of attributes: major-axis, minor-axis, center-of-mass-x, center-of-mass-y and angle-of-major-axis, as ellipsoids (blobs) for all regions from Figure 5 (left). Euclidean distance of the mean values for the attributes of all blobs is smallest to the single values of the indicated region (right).

4 Manual object analysis

For manual phase contrast image analysis ImageJ, a public domain, Java based image processing program was used (http://imagej.nih.gov/ij/docs/intro.html). Where necessary, the contrast and brightness of the images was adjusted.

For the low cell density images the cells bodies were manually outlined. By fitting an ellipse to the selected region, the angle between the primary axis of the ellipse and a line parallel to the X-axis of the image was computed, hence providing us with the single cells orientation (Figure 7). The average cell orientation was determined by averaging all the angles obtained from one image.

Figure 7

(Left) Manual analysis of single cell orientation using ImageJ. (Right) Application of the Fast Fourier Transformation analysis to the sample image. The red lines represent the average orientation of each 48 × 48 pixel square. The inaccuracies in the analysis, shown with the green arrows, could be explained by insufficient contrast in phase contrast image between the cells and the background.

The high cell density images were not possible to analyse using the manual delineation method since in such images it is very difficult to distinguish individual cell bodies. For that purposes the cell orientation was obtained by using an ImageJ macro, which was initially written for evaluation of F-actin orientation inside the cells. Shortly, the contrast of the images was enhanced and the cell orientation was extracted by 48 × 48 pixel sliding square analysis. Fast Fourier Transformation (FFT) performed in each square was followed with shape analysis of the FFT image by fitting the ellipse to the Fourier spectra and calculating the angle of the major axis of the fitted ellipse. Rotation by 90 degrees yielded the mean orientation of the analysed region. To visually ensure the accuracy of the data, obtained angle values were used to draw lines above the images (Figure 7). The average cell orientation was obtained by averaging the results obtained from individual analysed squares. The inaccuracies present in this analysis could be explained by insufficient contrast between the cells and the background in the phase contrast images. The obtained results of the regions with falsely computed orientation were disregarded in the quantitative analysis.

5 Discussion

The disadvantage of the manual image analysis is that it is subjective, hardly reproducible and extremely time consuming. It requires two separate methods to identify the cell orientation at low and high cell densities. In contrast to the explicit extraction of visually perceivable objects the analysis of uniformly distributed tiles yields incomprehensible artifacts and results since the orientations cannot be compared directly to an objective ground truth.

In the next step the automated and the manual approach will be quantitatively evaluated against an appropriate ground-truth with respect to accuracy (precision, recall) and coverage of single cell area segmentation and overall extraction rate.